FR2662524A1 - Process of three-dimensional visualisation by ray tracing with local interpolation - Google Patents

Abstract

The invention relates to all sectors which require means of three-dimensional visualisation and in which the data bases are described on acquisition as a stack of two-dimensional sections or directly in 3-D but with anisotropic sampling in the three dimensions. The process is therefore of relevance, among other things, to all medical imaging sources which provide this type of data. Conventionally, the isotropy of the bases described above is restored, before any manipulation or visualisation, through interpolation of sections. This process has the disadvantage of increasing the volume of data and hence of raising the size of the memory required in the processing systems. The invention is characterised in that no prior global interpolation is carried out on the anisotropic data base (7) before visualisation thereof by ray tracing. To do this, it follows the following procedure: - (8) tracing the ray in the anisotropic base, - (9) estimation of the position of the point of intersection of the ray with the subject to be visualised in the data base, - (10) interpolation in the vicinity of this estimated point, - (11) fine detection of the surface in this interpolated zone, - (12) shadowing as a function of the parameters contained in this zone.

Description

Three-dimensional ray tracing method with local interpolation.

 The invention relates to all domains which require three-dimensional visualization means (3-D) and in which the databases are described on acquisition as a stack of two-dimensional cuts or directly in 3-D but with anisotropic sampling in the three dimensions. The process therefore concerns, among other things, all medical imaging sources providing this type of data: It can be implemented on their host computer but also on any specialized or general purpose computer.

Currently three-dimensional visualization techniques can be identified by the type of representation of objects contained in the database they adopt. These are - surface representations (description of the surface of
the object by approximations of order 1 or more, or by
voxels - elementary volume usually cubic), - or voluminal representations, which contrary to the preceding
you keep all the information in the database
initial data.

 In each of these great frameworks, many methods have been proposed. The invention relates to the representation by volume and more particularly ray tracing techniques (as opposed to projection type techniques or octal shaft type or "octree").

Approaches based on ray tracing appeared in the 1970s in the synthesis of images: the search for special effects (multiple reflection, refraction, etc.), however, resulted in the complexity of their applications. Their use for rendering volumes only goes back three or four years. Visualization of volumes brings a major simplification, however, as realistic perceptions of databases will be favored over the search for special effects. Only then will the primary rays be calculated. These techniques are now perfectly tested and in the medical field we can refer to Drebin's articles (DREBIN R., CARPENTER L. and HANRAHAN P., Volume rendering,
ACM Com.Graph., 12, 4, Aug. 88, pp. 65-74), by Levoy (LEVOY M.,
Thesis, University of North
Carolina at Chapel Hill, 1989), by Robb (ROBB R. and BARRILLOT C.,
Interactive 3-D display and analysis, SPIE vol 939 Hybrid Image and
Signal Processing, 1988, pp. 173-202) and Trousset (TROUSSET.Y.,
Visualization of three-dimensional medical objects, Thesis, ENST Paris, September 1987 and French patent FR 2.613.509).

 All these works have brought improvements to the basic pattern of ray tracing, either in terms of the quality of representations, or in terms of computation time. The method, described later, has the advantage of reducing the memory capacities required for the 3-D representation without having to limit the other functionalities (geometrical transformations, cut by any planes, detection of sub-objects, qualities of the picture, ...).

In 3-D image visualization, a compromise must be found between three main factors: (a) the quality of the images and functions offered to the user, (b) the response times of the system knowing that interactivity
man-machine (even for large database sizes)
must be ensured, (c) the costs of the imaging stations.

 The proposed invention makes it possible to bring a significant gain on the third point (c) by decreasing in particular the size of the memory necessary for the processing systems.

 The peculiarity of most 3-D databases made from sections (or planes) 2-D is to present an axial anisotropy. That is, the resolution -the sampling in the plane is better than in the orthogonal direction.

 All the works, published to date, restore the isotropy of the volume by interpolating between cuts before any visualization or manipulation of the database. The total size of the database is thus multiplied by a factor N, with N # 2.

 It seems advantageous, to avoid this increase in the volume of data, not to practice a prior global interpolation and to directly view the content of the anisotropic base. A single publication (TROUSSET.Y., Visualization of three-dimensional medical objects, thesis, ENST Paris, Sept. 1987) proposes a scheme to economize such an increase but by constraining the problem so that the condition ( a) is no longer fully respected. Its method of viewing a three-dimensional database formed of stacked sections is characterized in that its incidence of vision is chosen so that the screen is placed perpendicular to the sections and that the rays follow a parallel geometry. Thus, each line of the screen corresponds to the image of the slice of a section. After a ray tracing on the anisotropic base, it is enough for it to interpolate lines between the lines of the screen. This method is generalizable neither in the case where the incidence of vision is any, nor in that of a conical geometry.

 The present invention aims to provide a ray tracing display method that does not require prior interpolation of the database and overcomes the aforementioned restrictions.

That is, - the initial anisotropic database remains unchanged, - only the visible parts (internal surfaces or volumes) to
the observer and regardless of the point of view chosen shall be
interpolated during visualization.

For this purpose, the invention proposes a method characterized by a sequence of operations to be performed for each radius and which is broken down into different phases (1) definition of the radius (trajectory) in the anisotropic base, (2) estimation of the position from the point of intersection of the radius with the
or the subjects to be visualized in the database, (3) interpolation on the neighborhood of this estimated point (the choice of
the interpolator is free, it only needs to be local), (4) search on the radius in the interpolated area of the point
intersection, (5) possible calculation of normal, (6) specification of image value (shading).

 The invention will be better understood on reading the description which follows, of a particular embodiment of the invention, given by way of non-limiting example. This example uses the basic principles of ray tracing and defines the essential notations. For the sake of simplification, the example is described in parallel geometry. Its principle, however, extends to future projections.

The description refers to the diagrams that accompany it in which
FIG. 1 is a diagram showing the two markers (screen
and objects) and the geometric relations between them,
FIG. 2 is a diagram showing the ray tracing in
parallel geometry,
FIG. 3 is a graph intended to illustrate the opera suite
ray tracing with local interpolation,
FIG. 4 is a diagram that poses the problem of intersec
between the department and the object described in the database,
FIG. 5 is a diagram illustrating one of the methods of
approximation of the position of the radius / object intersection,
FIG. 6 is a diagram illustrating a second method
approximate estimation of this position,
- Figure 7 shows histograms for evaluation
proposed methods.

A radius in an isotropic base is defined as follows
The 3-D scene has its own marker (1) (see Figure 1) called the Rob object mark. This depends on the arrangement of tomographic sections. Traditionally these are considered parallel to the plane (xobXyob) and stacked along the zob axis.

 The other coordinate system (2) -the screen mark or image mark tYc is linked to the visualization. The screen represents the plane (xec, yec). The scene is placed on the side of positive Zec.

 These two marks are linked by rotational and translational relations. These relationships depend on the vision parameters. The direction of the observation is defined by its spherical coordinates in the object reference (1) which form a rotation matrix [R] between the two marks. The origins of the landmarks are linked by a translation vector.

 When launching a ray, the relative position of a point of the ray relative to the object must be known continuously. This constantly requires to move from one benchmark to another. However, to visualize a database of voxels, it is more efficient to turn the screen around the object than to rotate the object itself.

It is therefore the screen / object transformations that will occur most often.

This amounts to applying the following formula
[Xob, Yob, Zob] = [R] [Xec, Yec, Zec] + T
The principle of ray tracing is as follows (see Figure 2). From the starting point (3) P0 (x0, y0, z0), the ray (4) is traversed from the screen (5) to the scene (6). At each increment, a screen / object mark transformation is performed. On the result of this one - that is to say the coordinates of the position of the radius in the reference object - are verified the conditions of intersection (by a function which we will call CONDITION ~ D ~ INTERSECTION). These can be based on density, local gradient, or other attributes of voxel. Once the intersection is found, the color of the pixel (xo, y0) is calculated.

The ray is a straight line in the screen mark. The parameters of this ray are
- its starting point Po (x0, y0, z0),
its director vector d [dx, dy, dz],
a variable t representing the position along the radius.

x(t) = x0 + doit y(t) = y0 + z(t) = z0 + dz.t
Si les rayons sont parallèles entre eux et à l'axe Zec t vaut [ 0,0,1 ] . His equation is written

x (t) = x0 + must y (t) = y0 + z (t) = z0 + dz.t
If the rays are parallel to each other and to the axis Zec t is [0,0,1].

x(t) = xo y(t) = y0 z(t) = z0 + t c'est à dire : [xec(t),yec(t),zec(t)] = [ x0,y0,z0+t ] où x0 et y0 sont les coordonnées du pixel de départ,
Les coordonnées du rayon dans le repère objet, pour un t donné, sont
[xob(t),yob(t),zob(t)] = [R]. [xec(t),yec(t),zec(t)] + T
= [R].[x0,y0,z0+t] + T
Cette opération se décompose en trois termes

xob(t) = x0.r11 + y0.r21 + (z0+t).r31 + tx yob(t) = x0.r12 + y0.r22 + (z0+t).r32 + ty
Zob(t) = x0.r13 + y0 r23 + (z0+t).r33 + tz
avec rjj les composante de la matrice [R] et ti les valeurs du vecteur T.The radius equation can be reduced to

x (t) = xo y (t) = y0 z (t) = z0 + t ie: [xec (t), yec (t), zec (t)] = [x0, y0, z0 + t] where x0 and y0 are the coordinates of the starting pixel,
The coordinates of the radius in the object reference, for a given t, are
[xob (t), yob (t), zob (t)] = [R]. [xec (t), yec (t), zec (t)] + T
= [R]. [X0, y0, z0 + t] + T
This operation breaks down into three terms

xob (t) = x0.r11 + y0.r21 + (z0 + t) .r31 + tx yob (t) = x0.r12 + y0.r22 + (z0 + t) .r32 + ty
Zob (t) = x0.r13 + y0 r23 + (z0 + t) .r33 + tz
with rjj the components of the matrix [R] and ti the values of the vector T.

 The transformation of the coordinates at each step of the ray slowed down the process. We will therefore increment the radius directly in the object reference.

Xob(ti) - xob(tj-1) = (ti-ti-1).r31
Yob(ti) - yob(ti-1) = (ti-ti-1).r32 zob(t;) - zob(ti-1) = (ti-ti-1).r33 Si nous considérons # #t = ti-ti-1 (pas du rayon), alors nous avons:
[xob(ti ) ,yob(ti) X zob(ti) ] = [ xo#(ti.i) ,yob(ti-1 ) , Zob(ti-1)
+ [ r31,r32,r33 ] . This increment is worth

Xob (ti) - xob (tj-1) = (ti-ti-1) .r31
Yob (ti) - yob (ti-1) = (ti-ti-1) .r32 zob (t) - zob (ti-1) = (ti-ti-1) .r33 If we consider # #t = ti-ti-1 (not radius), then we have:
[xob (ti), yob (ti) X zob (ti)] = [xo # (ti.i), yob (ti-1), Zob (ti-1)
+ [r31, r32, r33].

 As a reminder and with reference to the block diagram of FIG. 3 are the various operations of ray tracing on an anisotropic base (7) (the latter, in this example, will be composed of sections presenting a single anisotropy according to Zob): (8 ) launching the radius into the anisotropic base, (9) estimating the position of the point of intersection of the radius with the subject (s) to be displayed in the database, (10) interpolation on the vicinity of that estimated point, (11) ) fine detection of the surface in this interpolated zone, (12) shading according to the parameters contained in this zone.

The radius in the anisotropic base is defined as follows:
The sampling of the Rob landmark is no longer isotropic. If D represents the inter-slice distance and p is the size of a pixel of the original slices, the compression ratio (denoted e) of Rob along the zob axis is: e = D / p.

The isotropic resampling, by interpolation of the fictitious sections between two original sections, defines a new reference point called Riso. The transformation between R. and Rob is subjected to a scaling matrix E.

<Tb>

<SEP> 1 <SEP> 0 <SEP> 0 <SEP>
<tb> E <SEP> O <SEP><SEP> 1 <SEP> 0
<tb><SEP> 0 <SEP> Y <SEP> l / e <SEP>
<Tb>
The ray tracing, previously described, is done from the mark on on Rjso. The geometric transformations between the different marks are
[xiso, yiso, ziso] = [R].

[xob, yob, zob] = [E]. [XisosYisoXZiso]
The estimation of the position of the intersection of the radius with the subject (s) to be displayed in the database can proceed as follows:
See Figure 4; if for the original planes (13) the spokes (14) touch the object (15), for the "intermediate" planes a still non-existent surface must be detected.

 the interpolation on a volume neighborhood is greedy in computing time. The reduction of this neighborhood imposes an estimate of the position of the surface as accurate as possible.

Pa(x,y,Z) : point de détection approchée de la surface viso(xiso,yiso,ziso) : voxel de Rjso de position (xiso,yiso,ziso)
Rayon(xec,yec) : rayon lancé du pixel (xe#,yec) ESTIMATION(x,y,z) : module d'estimation de la position de la surface
Plusieurs modèles de la fonction ESTIMATION sont étudiés, qui font intervenir des critères géométriques et les propriétés physiques des structures.The detection module is as follows: the estimated intersection point is a voxel of the Riso coordinate system which belongs to the trajectory of the radius starting from (xec, yec) and which satisfies a criterion called
ESTIMATE. That is to say

Pa (x, y, Z): approximate detection point of the viso surface (xiso, yiso, ziso): position Rjso voxel (xiso, yiso, ziso)
Radius (xec, yec): ray of the ray (xe #, yec) ESTIMATION (x, y, z): modulus of the position of the surface
Several models of the ESTIMATION function are studied, which involve geometric criteria and the physical properties of the structures.

the model
This model, which refers to Figure 5, is based on the original discretizations of the database. A ray tracing on an anisotropic base is, in fact, a ray tracing (16) on an isotropic basis for which each trench (17) returns to its original thickness (18). Everything happens as if in the reference Riso the voxels of the intermediate planes took the value of the original planes. We can therefore detect the voxels (18) of the surface in a conventional manner that is to say as soon as a ray (16) penetrates (19).

The approximate detection condition is written

et la partie entière de

CONDITION ~ D ~ INTERSECTION is the function mentioned above that checks whether the voxel belongs to the object or not. The candidate voxel has for address

and the whole part of

This translation of 0.5 makes it possible to distribute the thickness on both sides of the original section.

- 2nd model
This model refers to Figure 6. The radius (20), calculated in real numbers, is surrounded by other rays formed by the integer coordinates of the initial radius. These rays, which we may call auxiliary rays, necessarily affect the known areas of the object. The coordinates of the approximate position of the intersection of the radius are determined by interpolation of the intersecting positions of the auxiliary radii with the object.

 Only the interpolation according to Zob is decisive. Auxiliary rays can be reduced to the number of two: the one above (z upper integer) and the one below (z lower integer).

 The estimation proceeds as follows: the intersections of the upper (21) and lower (22) auxiliary radii are searched for by the CONDITIONP 'INTERSECTION (x, b, y, b, z, b + 1) for the upper auxiliary radius and by CONDITION ~ OF INTERSECTION (xOb, yOb, zOb) for the lower auxiliary radius. The point of the estimated area is represented by the intersection (23) of the radius (20) and the line (24) defined by the two intersection points found (21) and (22).

- 36 models
The first model can be refined. The physical parameter (density, ...) associated with the detected point is calculated by linear interpolation of the values of the surrounding cuts. If this value meets the intersection conditions, this point is retained, otherwise the radius is incremented.

e : rapport de compression d : distance entre

coupes originales selon z

with

e: compression ratio d: distance between

original cuts according to z

The local interpolation module is as follows
All interpolation techniques from local data are usable whether they are of the first order (linear interpolation) or more elaborate (spline functions, ...). This list is not limiting, the only restriction is the local character of the interpolation.

The size and shape of the area surrounding the estimated intersection point and in which the interpolation will be conducted may be arbitrarily defined. However, different criteria allow an adjustment of these parameters (this list is not exhaustive) - the accuracy of the detection model, - the ratio e (if e is low -2 or 3- the estimated position is
close to the "real" surface), - the chosen interpolation technique. It turns out that the points
estimated by the detection operator are close to the surface of
the object after linear interpolation. If the interpola technique
localization is linear, the region can be summed up
neighbors (related) of the detected voxel. For ordering techniques
this region must be extended, - the direction of vision. The region must be extended according to the direction
Lion of islon (that of the ray), the method of calculation of the shading (of the neighbors who are there
required).

The fine detection of the surface and the shading take place thus
The different voxels, coming from the interpolation module, are stored in a small auxiliary base. The ray, transposed in this auxiliary and isotropic base, proceeds in a classical way (search for the surface, shading - as a function of the gradient of the surface, the distance from the screen, sources of light, the nature of the area, .. ...).

The validation of the method was done on two levels - the subjective appreciation of the differences between images by
visual comparison, - the quantitative analysis of the values resulting from the proposed procedures.

 The latter could be performed directly in the interpolated 3-D space after detection. The local value of the comparisons and the complexity of implementation have led to a preference for comparing the projected images from different points of view. The procedure applied is as follows: 1) acquisition of 3-D anisotropic databases; 2) global interpolation and visualization at a particular angle of incidence; 3) for the same point of view, calculation of the image by local interpolation on the ray; 4) point-by-point subtraction of the images obtained in (2) and (3) and statistics on the differences.

 This evaluation was made for variable interpolated section numbers, using the previously described detection methods and different interpolation schemes. The results after local reparcours of the radius do not have significant difference, only the approximate detection points were used, which with interest to simplify the method.

 For one of the databases tested, based on 128 X 128 X 136 voxels (isotropic volume) resolution, the removal of every other cut was performed. The database, rendered anisotropic, was visualized either after global (linear) interpolation or by ray tracing with local interpolation. Figure 7 shows the histograms of the differences between the original image (starting isotropic database) and the images produced after global interpolation (25) or local interpolation (1st detection model (26), 2nd (27) and 3rd (28)). It shows that the 3rd model of detection gives the best results and that these are extremely close to those of the global interpolation.

 Naturally, the present invention is not limited to the description that has just been made. In particular, it can be applied to all techniques for accelerating the ray tracing process. The advantage of this method is that it can be generalized or inserted in any sequence of operators along the radius operators local distribution of physical parameters or simply geometric criteria.

Claims (8)

claims  1 - Method for three-dimensional visualization of the content of a volume database having one or more axial anisotropies (different sampling along the three axes), according to any incidence, using the principle of ray tracing and characterized in that it comprises the following operations - definition of the radius (trajectory) in the anisotropic base, - estimation of the position of the point of intersection of the radius with the  or the subjects to be visualized in the database, - interpolation on the neighborhood of this estimated point (the choice of  the interpolation is free, it only needs to be local), - search on the radius in the interpolated zone of the point  intersection, - possible calculation of the normal at this point and specification of the value in the image plane (shading).  2 - Process according to claim 1, characterized in that the position of the point of intersection of the radius with the subject or subjects contained in the anisotropic base is estimated when the ray penetrates a voxel (anisotropic) of this base, voxel which responds to desired intersection conditions.  3 - Process according to claim 1, characterized in that the position of the point of intersection of the radius with the subject or subjects contained in the anisotropic base is estimated by a method which makes use of the following principle: the radius, calculated in real number , is surrounded by other rays formed by the integer coordinates of the initial ray, that these rays, which we may call auxiliary rays, necessarily affect the known zones of the object; thus the coordinates of the approximate position of. the intersection of the radius are determined by interpolation of the intersections positions of the auxiliary radii with the object.  4 - Process according to claim 3, characterized in that the position of the point of intersection of the radius with the subject or subjects contained in the anisotropic base is estimated by a method which makes use of the fact that if the database only presents a single axial anisotropy only the interpolation in this direction is decisive; thus the auxiliary rays can be reduced to the number of two: those which are on both sides of the ray in this direction; the estimation proceeds the principle according to the intersections of the two auxiliary rays are sought, the point of the estimated surface is represented by the intersection (23) of the radius (20) and of the line (24) defined by the two points d intersection found (21) and (22).  5 - Process according to claim 2, wherein the physical parameter associated with the point found according to the method of claim 2 is calculated by linear interpolation of the values of the surrounding cuts; the position of the point of intersection of the radius with the subject or subjects contained in the anisotropic base is estimated if the value of the physical parameter meets the intersection conditions.  6 - Process according to any one of claims 1 to 5, characterized in that the estimated position of the intersection between the radius and the subject or subjects contained in the database is considered as the exact position of this intersection.  7 - Process according to any one of claims 1 to 6, characterized in that the size and shape of the area surrounding the estimated intersection point and in which the interpolation will be conducted are defined either arbitrarily or in one or more criteria of the following list - the accuracy of the detection model, - the relationship (s) characterizing the anisotropy, - the chosen interpolation technique, - the direction of vision, - the shading calculation method (of the neighbors who are there  required).  8 - Process according to any one of claims 1 to 7, characterized in that the different voxels, from the interpolation module, are stored in a small temporary auxiliary base and the radius, transposed in this auxiliary base and isotropic, proceeds in a conventional manner to find the object or objects to be displayed and to specify the value of the image plane (shading).