WO2016030692A1

WO2016030692A1 - Method and apparatus for modelling non-rigid networks

Info

Publication number: WO2016030692A1
Application number: PCT/GB2015/052495
Authority: WO
Inventors: Serkan ÇIMEN; Alejandro Federico FRANGI
Original assignee: The University Of Sheffield
Priority date: 2014-08-29
Filing date: 2015-08-28
Publication date: 2016-03-03
Also published as: GB201415307D0

Abstract

Embodiments of the present invention provide a computer-implemented method of determining three-dimensionallocations of a non-rigid network, comprising receiving from a detector data indicative of a plurality two-dimensional (2D) projections through the network, identifying a plurality of network landmarks in the 2D projections, determining a correspondence between the network landmarks and landmark points of a three-dimensional (3D) model comprising a plurality of landmark points defining a surface associated with the network, and adapting, based on the correspondence, one or more parameters of the model.

Description

Method and Apparatus for Modelling Non-Rigid Networks

Background It is sometimes desired to investigate a non-rigid network, such as in medical applications where the non-rigid network is a biological network. For example when investigating coronary artery disease (CAD) it is useful to describe a coronary artery tree. A considerable amount of information is lost in 2D projections such as X-ray images of such non-rigid networks, e.g. owing to foreshortening and overprojection, and it is difficult to observe dynamic variations of the network.

It is an object of embodiments of the invention to at least mitigate one or more of the problems of the prior art. Summary of the Invention

According to aspects of the invention there is provided methods and apparatus as set forth in the appended claims. A method of determining four-dimensional locations of a non-rigid biological network, comprising receiving data indicative of a plurality of time stamped 2D projections through the network, identifying a plurality of network landmarks in two- dimensional (2D) projections, determining a correspondence between the network landmarks and landmark points of a four-dimensional (3D+t) model comprising a plurality of landmark points defining a dynamic surface associated with the network, and adapting, based on the correspondence, one or more parameters of the model.

An apparatus for determining three-dimensional locations of a non-rigid network, comprising a processor and a memory storing computer-executable instructions, wherein the instructions, when executed by the processor, are arranged to perform a method comprising the steps of receiving data indicative of a plurality 2D projections through the network, identifying a plurality of network landmarks in the two- dimensional (2D) projections, determining a correspondence between the network landmarks and landmark points of a three-dimensional (3D) model comprising a plurality of landmark points defining a surface associated with the network, and adapting, based on the correspondence, one or more parameters of the model.

Brief Description of the Drawings

Embodiments of the invention will now be described by way of example only, with reference to the accompanying figures, in which:

Figure 1 shows a method according to an embodiment of the invention;

Figure 2 shows an illustration of a plurality of images through an object according to an embodiment of the invention;

Figure 3 shows a body and a plurality of image planes intersecting the body according to an embodiment of the invention;

Figure 4 illustrates an image 400 of a network according to an embodiment of the invention; and Figure 5 illustrates a performance of an embodiment of the invention. Detailed Description of Embodiments of the Invention

Figure 1 illustrates a method 100 according to an embodiment of the invention. The method is a method of modelling a non-rigid network. In the described exemplary embodiment the non-rigid network is a coronary artery tree. It will be realised, however, that biological non-rigid networks of other types of may also be modelled, such as other biological networks (e.g. airway trees, cerebral vascular trees, neuronal networks, etc.). For example a network of airways surrounding a biological structure, such as the lungs may be modelled. Embodiments of the present invention are useful in modelling non-rigid networks, which are associated with an object or surface where the network may be attached to, follow or surround, the surface or object. For example the object may be the heart or its ventricular chamber where the non-rigid network is a coronary artery tree which at least partly surrounds the heart. Since the non-rigid network is associated with the object, the shape and movement of the object influences a shape of the network. Particularly where the object is heart which has distinct phases of movement, it can be appreciated that the coronary artery tree also moves in association with the heart.

Embodiments of the invention utilise a model which is representative of a shape or the surface or object. In embodiments of the invention the model is formed by a plurality of points which model the shape of the object. The model may model two or more factors of the object. In some embodiments the model accounts for inter-subject variation and temporal variation of the object. Inter-subject variation is variance in the shape of the object between instances i.e. persons of a class population, whilst temporal variation is a variation in the shape of the object over time i.e. changes in shape due to phase of a cardiac cycle or any other source of motion. The model is used to constrain shape and motion of the non-rigid network being modelled.

Referring to Figure 1, the method comprises a step 1 10 of providing a model of an object associated with the network. In some embodiments of the invention a model of a surface of the object upon which the network lies (e.g. the ventricular epicardium, in this example) is provided, although it will be realised that other structures, particularly biological structures, may be modelled. The model is a bilinear model which is suitable to model variations due to two independent factors. The two independent factors may be referred to as style (shape from patient or subject) and content (time or phase). In other embodiments of the invention, a larger number of factors might be considered by using multi-linear object decompositions instead of bilinear object decompositions. The essence of any object decomposition in any embodiment is that it helps to regularise the dynamic reconstruction of the network by incorporating prior knowledge associated with the three-dimensional object and its dynamics.

As shown in Figure 2, the model 200 comprises a plurality (N) of landmark points 210, 220, 230, 240 defining locations on a surface of the modelled object. It will be realised that the number of points 210, 220, 230, 240 illustrated in Figure 2 is not restrictive. The points may be joined to form a mesh indicative of the shape of the object, such as the heart ventricles or the lungs. Each point 210, 220, 230, 240 of the model is associated with a probability value p_n. The probability value is indicative of a probability of the respective point of the model corresponding to a location in the non-rigid network. For example, where the model 200 is indicative of a heart the probability value may be determined with respect to an atlas of the heart. The atlas of the heart defines exemplary locations for the network, such as coronary arteries with respect to the heart. The probability value is indicative of the point corresponding to an artery point in one embodiment.

In one embodiment the probability value p_n is based upon a distance of the point 210, 220, 230, 240 to a centreline of a coronary artery. The probability may be calculated as:

- ²

exp

Pn

(ψ)²

∑n=l ^exP ^ζ

Where N is a number of points 210, 220, 230, 240, d„ is a distance of a respective point n to a coronary artery on a surface of the heart and ζ is a constant which avoids strict correspondences due to p_n values and provides flexibility to compensate for anatomical variability. Thus it can be appreciated that the probability value p_n provides prior information to associate the network with the model 200.

The model 200 may be constructed based upon a training set of temporally and spatially aligned cardiac ventricle surfaces of S subjects (patients) in C cardiac phases. Each training ventricle surface is represented by N landmark points in a <i-dimensional Euclidean space, which in some embodiments of the invention is 3 -dimensional space (d=3). The landmark points are concatenated to form a K(= N X d)dimensional observation vector, y^sc. By using the bilinear model 200, each element of y^sc is written as yj = a^sTW_kb^c where a^sand b^c denote bilinear model parameters of subject and phase which are / and J dimensional vectors, respectively. In one embodiment a^sis a vector indicative of parameters defining the shape of the object based on the subject and b^c is a vector of parameters defining the shape of the object based on phase. The vector a^smay comprise more parameters than the vector b^c . In one embodiment the vector of parameters a^s comprises 99 parameters and the vector of parameters b^c may comprise 7 parameters, although it will be realised that other numbers of parameters may be chosen. W_k is an / X / matrix determining an interaction of two factors. One way to select the number of parameters is by selecting the number that ensures the model explains a given percentage (e.g. 99%) of the variance in the population. The model 200 is constructed by a training process which determines the parameters a^sand b^c for each subject class s and phase class c and the interaction matrices W_k . The parameters may be determined by minimising a total squared error. An iterative method may be used, such as based upon singular value decomposition (SVD). The model training results in an / X 5 matrix A = [a¹ ... a^s], J X C matrix B = [b¹ ... b^c] and IK X / matrix W

It can be appreciated from above that in embodiments of the invention a four- dimensional model (three dimensions + time (t)) is used to model an obj ect associated with the non-rigid network. In some embodiments the model parameters have been determined based upon a training set of object measurements of real-world objects. The model 200 may be stored in a memory of an apparatus performing the method 100 according to an embodiment of the invention. In some embodiments of the invention N is at least 1000 and may be around2000 landmark points, although it will be realised that other values of N may be chosen appropriately depending on the type and complexity of the network. The model may comprise a plurality of phases indicative of the shape of the object during each phase. Each phase may be associated with a numeric identifier indicative of the phase.

In step 120 data indicative of a plurality 2D projections through the network is received. The data may be received over a computer network. The data may be received from a detector apparatus arranged to determine the plurality of projections. In some embodiments the apparatus may be an apparatus for X-ray rotational angiography. The apparatus for rotational angiography comprises a portion known as an arm, or more specifically a C-arm, which rotates around a patient and outputs the plurality of 2D projections in the form of a series of images at each of a plurality of rotation angles. A radiation source emits radiation at each rotation angle and a detector receives radiation transmitted through the body 300 to determine the image. A projection matrix Pf is associated with each image which provides a relationship between three dimensional points and the two dimensional images recorded by the detector, as will be appreciated.

Referring to Figure 3 there is illustrated an object 300 which may be a patient's body comprising a heart. The X-ray rotational angiograph rotates around the body 300 to determine a plurality of images intersecting the body 300 at each of a plurality of angles 320, 330, 340. It will be realised that the number and respective angle of each of the images is not limited by those illustrated in Figure 3. The images also intersect the network of interest, such as the coronary artery tree. There are F images through the network where /is indicative of a respective image.

As explained above, where the network is associated with a moving object, that is an object which may be stationary in position but variable in shape and/or size, such as the heart, over a period of time during which the series of images are produced, the object may change shape and/or size. In particular whilst the arm is rotating around the body 300 the heart may repeatedly move through various cardiac phases.

In embodiments of the invention the images are each associated with a respective phase of movement of the object. In one embodiment one or more measurements are recorded whilst determining the images to assign each image a timestamp representative of the respective dynamic phase of the object. For example each image may be assigned to a respective cardiac phase. The timestamp may be derived from measurements associated with the object, such as of electrocardiogram (ECG) or blood-pressure (BP) signal. For example, the ECG is indicative of electrical activity of the heart. Based on one or more features of the ECG signal, such as RR intervals (where the interval is between consecutive R waves, as will be appreciated) of the ECG signal, each image is assigned to a respective cardiac phase. The image may be assigned to the cardiac phase by assigning a numeric identifier to the image indicative of the phase. The identifier may be between first and second numbers such as between 0 and 1 although it will be realised that other numbers may be chosen. An interval of the numbers may the same as a range of the numeric identifiers associated with phases of the model, such that the identifier associated with the image may be used to identify a corresponding phase of the model 200. In step 130 a plurality of network landmarks in the 2D projections are identified. The network landmarks correspond to characteristic locations in the network identifiable in each of the images determined in step 120. In some embodiments between 100 and 200 network landmarks may be identified in each image, although it will be realised that this is merely exemplary and that other numbers of network landmarks may be used particularly depending upon a size of the network. The number of network landmarks may be M and, in some embodiments M«N. In some embodiments may be at least 50, at least 100 or between 100 and 200 network landmarks in each image, although it will be realised that other numbers of network landmarks may be identified. Identification of the network landmarks may be performed by a variety of image analysis methods including feature extraction and interest point detectors, as will be appreciated.

Figure 4 illustrates an exemplary image 400 of a network comprising a plurality of visible branches, such as coronary arteries. The image may have been determined using the X-ray rotational angiography apparatus described above. As a result of step 130 a plurality of network landmarks 410, 420, 430, 440, 450 are identified at locations in the image corresponding to points within the network. It will be appreciated that the network landmarks do not have to be at junctions in the network and may be at intermediate locations along branches of the network.

In step 140 a correspondence between the network landmarks 410, 420, 430, 440, 450 determined in step 130 and the model 200 is determined. The correspondence is determined between the two dimensional network landmarks 410, 420, 430, 440, 450 identified in the image and projections of the landmark points 210, 220, 230, 240 of the model. Because the landmark points 210, 220, 230, 240 of the model are in three dimensions defining the shape of the object, they are projected to a plane of the image in two dimensions. In embodiments of the invention the network landmarks 410, 420, 430, 440, 450 are assumed to be of fixed position, whilst parameters of the model 200 are altered in order to vary a location of the projections of the landmark points 210, 220, 230, 240 in the plane of the image. The correspondence defines a transform between two point sets, namely the network landmarks 410, 420, 430, 440, 450 and the projections of the landmark points 210, 220, 230, 240 in the plane of the image. Based on the correspondence parameters of the model are determined, as will be explained.

The correspondence is represented in some embodiments as a matrix G. The correspondence matrix G may by M x N in size. The correspondence matrix may indicate, for each landmark point 210, 220, 230, 240 a likelihood of corresponding to each network landmark 410, 420, 430, 440, 450. As can be appreciated from Figure 1, steps 140-180 are iteratively performed. That is, steps 140-180 are repeatedly performed until a predetermined condition is met. Steps 140-180 determine a registration between the two point sets. During iterations of the method 100 the position of landmark points 210, 220, 230, 240 in the model 200 associated with the network landmarks 410, 420, 430, 440, 450 is gradually refined.

In embodiments of the invention a determination of registration between the two point sets is determined based on a Gaussian mixture model. In embodiments of the invention, it is assumed that projections of landmark points 210, 220, 230, 240 in the model 200 specify Gaussian cluster centres whilst the network landmarks 410, 420, 430, 440, 450 describe a spatial distribution of points in two dimensions. A mean of each Gaussian cluster is given by a projected location of a corresponding landmark point 210, 220, 230, 240 in the model 200. Thus the mean of each Gaussian distribution is determined by parameters of the model 200. Each Gaussian distribution is further defined by a standard deviation in each of the two dimensions (x, y).

As noted above, each landmark point 210, 220, 230, 240 in the model 200 is associated with a probability value. The probability values may be used as a weight for each Gaussian distribution. In this sense the weight may be referred to as a mixing coefficient, where a sum of all weights is unity. The Gaussian mixture model is defined as a sum of Gaussian distributions multiplied by the corresponding weights. Steps 140-180 aim to minimise a difference between the fixed distribution of network landmarks 410, 420, 430, 440, 450 and a mixture (or weighted sum) of Gaussian distributions.

A log-posterior energy function is used in some embodiments to find parameters of the model 200 and the correspondence matrix G. The parameters of the model found are a and B which are a vector and matrix indicative of subject and phase, respectively. It will be noted that B is a matrix due to storing model parameters for multiple cardiac phases i.e. a vector b^c for each of a plurality of cardiac phases. Furthermore, in some embodiments, a rotation matrix R and translation vector t are also determined which are indicative of rotation and translation between the model 200 and a coordinate system of the apparatus used to produce the images, such as an X-ray coordinate system. The energy function may be a weighted sum of a distance term Ea_st and a regularisation term E_reg. The distance term measures a sum of squared distance between 2D observations, namely the network landmarks 410, 420, 430, 440, 450, and 2D projections of points 210, 220, 230, 240 of the model 200.

The distance term may be:

where x_m\s the m observation or network landmark from the/ image and x„is a 2D f

projection of the model point estimate, , which may be given as: [W^VTa ^VTb + i

1

Where is the projection matrix which may be extracted from X-ray image tags and superscript VT is indicative of a vector transpose operation. W is a 31 x J matrix formed by rows of W which corresponds to one of the landmark points, b is a column of B corresponding to a phase of the ^h image. Where the phase of an image does not coincide with a discrete phase of the model training set, embodiments of the invention use estimates at neighbouring discrete time points and interpolate to find .

The regularisation term E_reg may be based upon prior model parameters. In one embodiment the regularisation term is defined as a negative log-likelihood of the prior distribution of model parameters. The model parameters learned during training of the model may be used to perform a kernel density estimation of components of the parameter vector a^sand matrix B. In one embodiment the regularisation term may be defined as:

Where h denotes a kernel width and /c(-) is a Gaussian kernel. The kernel width acts as a smoothing parameter for the kernel. In the case of a Gaussian kernel h may be interpreted as a standard deviation of a Gaussian probability density distribution.

Total energy is therefore defined as E_tot = E_dist + xXE_reg where λ is a weighting between the two energy terms and τ is an annealing parameter. In some embodiments annealing parameters ¾ Tf_m and x_up are used, where ¾ is an initial annealing parameter at a start of the method, Tf_m is a final annealing parameter which may be used to define an end of the method, and x_up is used to reduce the initial annealing parameter with increasing iterations. The annealing parameters ¾ Tf_m and x_up may be 500, 3 and 0.95, respectively, although it will be realised that these are merely exemplary.

The annealing parameter τ may be reduced in increasing iterations of steps 140-180 to reduce a search space for the method. The search space may be reduced by controlling the Gaussian distributions, such as by reducing a standard deviation of the Gaussian distributions in the Gaussian mixture model. The standard deviation may be reduced isotropically i.e. the same standard deviation value in x and y directions or anisotropically i.e. using different standard deviations in the x and y directions. As the standard deviation is reduced only network landmarks closer to the landmark points are assigned a higher correspondence.

In step 140 the correspondence matrix G is updated. In one embodiment the correspondence matrix is updated by:

Where p_n is the mixing coefficient, discussed above, above for the nth Gaussian cluster, τ is the annealing parameter and c is a constant which, in one embodiment, is given by: w \/2πτ

C^~ 1 - w M

Where w is weight for a uniform distribution. A uniform distribution is generally added to the Gaussian mixture model in order to model outlier points. It is a good practice to determine the effect of the uniform distribution depending on the expected outlier percentage in the data. This effect is controlled by the weight parameter w. The uniform distribution is weighted by w and all the other Gaussians distributions are weighted by (l-w) to determine that the result is a valid probability.

Steps 150 and 160 represent a process of discarding outlier data points. In step 150 it is determined whether G_mn is less than a threshold β. In some embodiments β is set to 0.05 although it will be realised that other values may be used. If the G_mn is less than β then ^and its corresponding Gaussian distribution is discarded in step 160. Thus a number of Gaussian distributions is gradually reduced by steps 150, 160 in iterations of the method.

In step 170 parameters of the model 200 are updated. The parameters of the model are updated by minimising the total energy E_tot as discussed above. In one embodiment of step 170 two minimisations of the total energy E_tot are performed. One minimisation is for the parameters a, R, t indicative of the subject, rotation and translation. Another minimisation is for the parameter matrix B. Advantageously performing separate minimisations ensures independence of subject and phase.

In step 180 it is determined whether to terminate further iterations of steps 140-180. The determination of whether to terminate the method 100 may be based upon one or more criteria, for example whether a predetermined number of iterations have been completed. In one embodiment iterations of steps 140-180 are terminated based upon the annealing parameter τ and a predetermined threshold. Firstly in step 180 the current value of the annealing parameter is reduced by a predetermined value, for example by multiplying τ by x_up. In one embodiment if the annealing parameter τ is less than a threshold value Tf_m then the method ends. If, however, the annealing parameter τ is not less than the threshold then it is reduced by a predetermined value and steps 140-190.

In order to verify embodiments of the invention a model comprising training surface meshes describing the left ventricular epicardium (N = 2044) were obtained using an atlas based segmentation algorithm from 134 retrospectively ECG-gated multi-slice CT images. These meshes were temporally aligned to compensate for heart rate differences between patients.

In order to quantitatively evaluate embodiments of the invention, synthetic X-ray rotational angiography data was generated using a single left coronary artery geometry from a 4D XCAT phantom. The cardiac cycle was set to be 1000ms, with no respiratory motion present. X-ray imaging parameters, including projection matrices, number of images (117) and frame rate (30 fps) were derived from a clinical dataset. A total of 208 corresponding coronary artery points were generated and projected to generate the 2D observation points in 117 projection images.

The phantom data was used to evaluate a robustness of the method 200 to missing data and noise. In a first experiment (Experiment 1) 0% to 50% of 2D observation points which describe the projected 2D trajectory of each tracked coronary artery point along the sequence of projection images were randomly removed. In a second experiment (Experiment 2) 2D observation points were removed according to a sampling scheme which mimics the sparsity associated to missing detail in the detection of centerline points in the 2D views. To this end, a set of points, Pi, were selected which consisted of a starting point, bifurcations and end points. Similarly, a set of points was defined P₂ which are midpoints of the segments described by the points in the set Pi and defined a set P₃ as the set of point which are the midpoints of the segments described by the points in the set Pi and P₂. An embodiment of the method 200 was applied, given the points in the sets Pi (9: 13% of all points), PI P2 all points (16:35% of all points) and PI P2 P₃ (28:37% of all points). Performance of the method 200 was evaluated under uncertain measurements by adding zero mean Gaussian noise (σ = 0.25 to 1.25 mm) to the 2D points (Experiment 3). In all of the experiments, reconstruction errors were measured as root-mean- square errors in 3D between the reconstructed points and the true 3D positions.

The parameters of the algorithm are set empirically. The annealing parameters, X Tf_m and x_up were set to 500, 3 and 0.95, respectively. The threshold for rejecting a projected bilinear point, β, was set to 0.05 and the weighting of the regularization term λ was set to 2.5 x 10⁵. Finally, the weight for the uniform outlier cluster in the GMM, w, was set to 1.0 x 10^"8 since it was assumed that there were no outliers in the 2D observations.

Figure 5 illustrates results of these experiments. Figure 5(a) shows the results of the Experiment 1, Figure 5(b) shows the results of the Experiment 3, Figure 5(c) shows the results of the Experiment 1 for all cardiac phases, Figure 5(d) shows results of the Experiment 2 for all cardiac phases , Figure 5(e) shows qualitative results of the Experiment 1 for end-diastolic and end-systolic cardiac phases, and Figure 5(f) shows qualitative results of the Experiment 2 for Pi and Ρ^Ρ₂^Ρ₃. For the qualitative results, a ground truth centerline 510 of the coronary artery tree is shown in sold line. Reconstructed points are shown as spheres if a 2D observation is available in the corresponding image, and as triangles if the 2D observation is removed.

The results show that the 3D+t reconstruction performance of the method 200 stays stable even under 1.25 mm 2D observation noise, which is approximately 6-7 times of the pixel resolution of the rotational angiography (Figure 5(b)). Although there are some outliers, the results indicate that the method 200 is able to handle missing data (Figure 5(a)). In particular, the qualitative results show that the 3D reconstruction of missing data points are recovered satisfactorily. The accuracy of the algorithm quickly increases with addition of a small number of points to the base set, P, and continue to increase as more 2D observations points are added (Figure 5(c)).

It will be appreciated that embodiments of the present invention provide a method of determining three-dimensional locations of a non-rigid network. Advantageously knowledge about the location of the network improves processes associated with the network. For example, where the network is a biological network, such as a coronary artery network although not limited thereto, knowledge about the location of the network may improve surgery or other treatment. In some embodiments, knowledge about the location of the network may be improved by determining one or more images comprising an indication of at least a portion of the network. Thus the method may comprise steps of outputting image data indicative of an image of the non-rigid network based on the model. The method may further comprise a step of generating an image indicative of the non-rigid network based on the image data.

It will be appreciated that embodiments of the present invention can be realised in the form of hardware, software or a combination of hardware and software. Any such software may be stored in the form of volatile or non-volatile storage such as, for example, a storage device like a ROM, whether erasable or rewritable or not, or in the form of memory such as, for example, RAM, memory chips, device or integrated circuits or on an optically or magnetically readable medium such as, for example, a CD, DVD, magnetic disk or magnetic tape. It will be appreciated that the storage devices and storage media are embodiments of machine-readable storage that are suitable for storing a program or programs that, when executed, implement embodiments of the present invention. Accordingly, embodiments provide a program comprising code for implementing a system or method as claimed in any preceding claim and a machine readable storage storing such a program. Still further, embodiments of the present invention may be conveyed electronically via any medium such as a communication signal carried over a wired or wireless connection and embodiments suitably encompass the same. All of the features disclosed in this specification (including any accompanying claims, abstract and drawings), and/or all of the steps of any method or process so disclosed, may be combined in any combination, except combinations where at least some of such features and/or steps are mutually exclusive.

Each feature disclosed in this specification (including any accompanying claims, abstract and drawings), may be replaced by alternative features serving the same, equivalent or similar purpose, unless expressly stated otherwise. Thus, unless expressly stated otherwise, each feature disclosed is one example only of a generic series of equivalent or similar features.

The invention is not restricted to the details of any foregoing embodiments. The invention extends to any novel one, or any novel combination, of the features disclosed in this specification (including any accompanying claims, abstract and drawings), or to any novel one, or any novel combination, of the steps of any method or process so disclosed. The claims should not be construed to cover merely the foregoing embodiments, but also any embodiments which fall within the scope of the claims.

Claims

A computer-implemented method of determining three-dimensional locations of a non-rigid network, comprising: receiving from a detector data indicative of a plurality two-dimensional (2D) projections through the network; identifying a plurality of network landmarks in the 2D projections; determining a correspondence between the network landmarks and landmark points of a three-dimensional (3D) model comprising a plurality of landmark points defining a surface associated with the network; and adapting, based on the correspondence, one or more parameters of the model.

The method of claim 1, wherein the model is a spatio-temporal model of the surface and each of the 2D projections corresponds to a respective point in time.

The method of any preceding claim, wherein determining the correspondence comprises determining a projection of each landmark point of the 3D model.

The method of claim 3, wherein each landmark point is projected onto a plane of the projections through the network.

The method of claim 3 or 4, wherein the correspondence is based upon a Gaussian distribution associated with the projection of each landmark point.

The method of claim 5, wherein the projection of each landmark point defines a mean of the Gaussian distribution.

The method of any preceding claim, wherein the correspondence is a matrix indicating, for each landmark point, a likelihood of corresponding to each network landmark.

The method of any preceding claim, wherein the steps of determining the correspondence and adapting the one or more parameters of the model are iteratively performed.

The method of claim 8, wherein, in at least some iterations, one or more outlying landmark points are discarded.

The method of claim 8 or 9, when dependent on claim 5, 6 or any claim dependent thereon, a search space is successively reduced in iterations of the method; optionally the search space is reduced by reducing a standard deviation of at least some Gaussian distributions in at least some iterations of the method.

The method of any preceding claim, wherein the one or more parameters of the model are adapted based on an energy function.

The method of claim 11, wherein the one or more parameters are adapted to reduce energy associated with the correspondence.

The method of claim 11 or 12, wherein the energy function comprises a distance term indicative of a distance between the network landmarks and a projection of each landmark point.

The method of claim 11, 12 or 13, wherein the energy function comprises a regularisation term based upon the one or more parameters of the model.

The method of any preceding claim, comprising receiving data indicative of one or more measurements associated with the network.

The method of claim 15, comprising assigning an identifier indicative of a phase of movement of the surface to each projection based on the data.

17. The method of any preceding claim, wherein the model is a bilinear or a multilinear model.

18. The method of claim 13, wherein the model is associated with one or more parameters defining a shape of the surface according to a subject, and one or more parameters defining a shape of the surface according to phase.

19. The method of claim 14, wherein the model is associated with one or both of a vector a^s of parameters defining the shape of the object based on the subject, and a vector b^c of parameters defining the shape of the object based on phase.

20. The method of any preceding claim, comprising determining a parameter indicative of one or both of rotation and translation between the model and the projections through the network.

21. The method of any preceding claim, wherein each landmark point is associated with a probability value indicative of the probability of the landmark point corresponding to a location in the network.

22. The method of claim 21 when dependent upon claim 5 or 6 or any claim dependent thereon, wherein the probability value is used as a mixing coefficient associated with a corresponding Gaussian distribution.

23. The method of any preceding claim, wherein the model is a model of a human or biological network structure.

24. The method of claim 23, wherein the network is a network of vessels, arteries, veins or airways; optionally the network is a network of coronary vessels .

25. The method of any preceding claim, wherein the plurality of 2D projections are determined based upon medical or biological imaging.

26. The method of any preceding claim wherein the plurality of 2D projection images are provided from a X-ray rotational angiograph.

27. The method of any preceding claim, comprising outputting image data indicative of an image of the non-rigid network based on the model.

28. The method of claim 27, comprising generating an image indicative of the non-rigid network based on the image data.

29. Computer software which, when executed by a computer, is arranged to perform a method according to any preceding claim.

30. The computer software of claim 29 stored on a computer-readable medium.

31. An apparatus arranged to implement a method according to any preceding claim.

32. An apparatus for determining three-dimensional locations of a non-rigid network, comprising: a processor; and a memory storing computer-executable instructions, wherein the instructions, when executed by the processor, are arranged to perform a method comprising the steps of: receiving from a detector data indicative of a plurality two-dimensional (2D) projections through the network; identifying a plurality of network landmarks in the 2D projections; determining a correspondence between the network landmarks and landmark points of a three-dimensional (3D) model comprising a plurality of landmark points defining a surface associated with the network; and adapting, based on the correspondence, one or more parameters of the model.