CN101174331B - Maximum density projection generating method for medical image - Google Patents
Maximum density projection generating method for medical image Download PDFInfo
- Publication number
- CN101174331B CN101174331B CN2006100634318A CN200610063431A CN101174331B CN 101174331 B CN101174331 B CN 101174331B CN 2006100634318 A CN2006100634318 A CN 2006100634318A CN 200610063431 A CN200610063431 A CN 200610063431A CN 101174331 B CN101174331 B CN 101174331B
- Authority
- CN
- China
- Prior art keywords
- image
- mip
- projection
- steps
- picture frame
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Abstract
The invention relates to a calculated medical image maximum density projection generation method, which comprises the following steps: step A, for an initial image data volume formed by a plurality of parallel image frames, the projection of one of the image frames on a projection plane is used as the mapping template; step B, after translation processing, the projection of other image frames on the projection plane is overlapped with the mapping template, and medical image maximum density projection processing is performed; step A comprises the following steps: step A1, the planar parametersand MIP image coordinate parameters of the projection plane are created, and the projection of the image frames on the projection plane is used as the mapping template; step A2, a mapping table is created; step B comprises the following steps: step B1, medical image maximum density MIP processing to each point of the first image frame is performed according to the mapping table; step B2, for the subsequent image frames, a new temporary mapping table is formed; step B3, corresponding medical image maximum density MIP processing to the subsequent image frames is performed according to the new mapping table, and the MIP image is dynamically refreshed to generate the final MIP image. The invention has the advantages of small operation quantity and small middle storage space requirement, thereby meeting real-time requests.
Description
Technical field
The present invention relates to calculate, relate in particular to a kind of maximum density projection generating method for medical image.
Background technology
Maximum intensity projection technology MIP is as a kind of visual 3-D technology, because the advantage of its simple possible, be effectively applied to the visual work of the blood vessel, bone and these structures of soft tissue that from magnetic resonance imaging MR and computer tomography CT data, obtain always, on the diagnostic work of diseases such as the cerebrovascular, cardiovascular and tumor tissues, can play important booster action.
As shown in Figure 1, the essence of MIP algorithm is exactly to pass a volume M with a branch of parallel ray S from certain angle, finds out the maximal value on every ray, and the image of being made up of these maximal values is called the MIP image.
The MIP algorithm is the bigger algorithm of a kind of calculated amount, and its ultimate principle requires all pixel values on each bar projection ray all to be obtained and compare.When light from arbitrarily angled projection the time, because different with its storage order to the addressing order of the data point in the polylith picture frame, the actual calculation amount is very big.Therefore, speed and picture quality become the focal issue of MIP research.
Draw in order to reach high-quality real-time, interactive, adopt parallel processor sometimes, and the algorithm of picture quality and render speed compromise.Suitably reducing under the situation of picture quality, can reach real-time interactive and draw.Some scholars have proposed to improve render speed and the method that do not influence picture quality, these methods adopt spatial data structure that sorted data fields is encoded, in computation process, only handle the voxel of transparency in certain scope, and need not handle the voxel of " transparent fully " and " highly opaque ", can reduce operand.These data structures comprise Octree (Octrees), pyramid (Pyramids) and k2d tree (k2dtrees) etc., these methods or pretreatment time are long, need very big redundant space, when direction of observation changes, then need again spatial data field to be done single treatment, cause to wait as long for.
In the prior art, adopt a kind of wrong cutting to separate algorithm (Shear-warp) sometimes, the Shear-warp algorithm is to decompose along the ray cast direction, to simplify the projection from the object scanning space to image space.Data fields from object space coordinate system transformation to a middle coordinate system, is called mistake and cuts body space (Sheared object space).In mistake was cut object space, all direction of visual lines were all cut into slices perpendicular to the voxel layer.As shown in Figure 2, object space is cut in expression from the object space to the mistake parallel projection conversion, horizontal line is represented the dropping cut slice of data volume, these sections are all perpendicular to direction of visual lines after the conversion.This conversion is translation transformation (as along x and y direction) in fact.In mistake is cut object space, voxel layer section (Voxels slices) can be projected to a middle transition image, again the middle transition image transformation be obtained final image.
As shown in Figure 3 and Figure 4, as follows based on the MIP image processing method of Shear-warp algorithm:
(1) determines primary axis according to direction of visual lines, data volume is transformed to mistake from data space cut body space, and, cut the section of Shear transfer pair by mistake and resample according to wrong contact transformation matrix.
(2) in wrong tangent space, the synthetic section that resamples.
(3) mistake is cut voxel and project to intermediate image plane.
(4), intermediate image is transformed to final image by decomposing the Warp conversion according to the deformation transformation matrix.
This mistake is cut the projection in the body space, and the following synthetic characteristic of section that is used to simplify is arranged:
Characteristic 1: the sweep trace of intermediate image is parallel with the sweep trace of each wrong slice.
Characteristic 2: the voxel of same section has identical scaling.
Characteristic 3: for parallel projection, all sections all have identical scaling, if elect scaling as 1, then the pixel of the sweep trace of Qie Pian sweep trace and intermediate image is an one-to-one relationship.
In the projection of section is synthetic, utilize character 1 can handle the sweep trace in wrong tangent space and intermediate image space simultaneously.Section is synthetic adopts back projection's method or sciagraphy forward.
Also must could form final image by the synthetic intermediate image that obtains of section, generally adopt the affined transformation of two-wire shape filtering to obtain final image through distortion.
Though existing Shear-warp algorithm improves on computing velocity, still has following shortcoming:
One, need carry out pre-service, promptly generate raw data field, normal vector field and gradient fields according to data volume, the processing time of these operations often surpasses the drafting time.
Two, the synthetic and intermediate image of section need take storage space, because the shear-warp algorithm can not guarantee that the intermediate image of current generation has maximal density in the synthetic process of section, could determine after must traveling through all sections by the time, produce to wait for end product so need to keep intermediate image.
Therefore, existing Shear-warp algorithm operation quantity is big, and required storage space is big.
Summary of the invention
The object of the present invention is to provide a kind of high efficiency maximum density projection generating method for medical image, to solve operand and the big problem of required storage space in the prior art.
Maximum density projection generating method for medical image of the present invention is:
A, for the raw image data body that constitutes by polylith parallel image frame, with the projection of picture frame on the projecting plane wherein as the mapping template;
B, the projection of other picture frame on the projecting plane are benchmark with described mapping template, handle through translation, and be overlapping with described mapping template, carries out the medical image maximum intensity projection and handle.
Described steps A comprises the steps:
A1, set up the plane parameter on projecting plane, set up corresponding M IP image coordinate parameter, with the 1st projection of picture frame on the projecting plane as the mapping template;
A2, generation mapping table, the pairing two-dimentional MIP image coordinate location of three-dimensional space position of data volume in the described mapping table reflection picture frame.
Described steps A 1 comprises the steps:
A11, determine the parameter on projecting plane, on the projecting plane, choose x direction of principal axis and the y direction of principal axis of two perpendicular directions as the MIP image;
A12, get the longest straight line in the raw image data body as the length of side of MIP image;
A13, the center of getting the intersection point conduct final MIP image of central point on the projecting plane of raw image data body.
Described steps A 2 comprises the steps:
Three dimensional space coordinate behind the spot projection in A21, the calculating chart picture frame;
A22, the volume coordinate after projection conversion is the two-dimensional coordinate in the MIP image;
A23, successively calculate in the 1st picture frame the two-dimensional coordinate in the somewhat pairing MIP image, and be stored in the mapping table.
In the described steps A 22, in the conversion of the two-dimensional coordinate in the MIP image, insert correction coefficient K
0
Described step B comprises the steps:
B1, according to described mapping table, each point of the 1st picture frame is made medical image maximal density MIP handles;
B2, for follow-up picture frame, with the mapping template be benchmark, each value in the described mapping table is made migration processing, form interim new mapping table;
B3, the new mapping table of basis are made corresponding medical image maximal density MIP to follow-up picture frame and are handled, and to MIP image dynamic refresh, produce final MIP image.
Among the described step B3, produce before the final MIP image, described MIP image is carried out the interpolation operation.
Gridiron pattern translation interpolating method is adopted in described interpolation operation.
Described MIP image is further carried out the operation of Gauss's smothing filtering.
Also comprise the steps A0 before the described steps A:
Axial anastomose property between A0, judgement projecting direction and the raw image data body, carry out following operation:
A01, if axially coincideing of projecting direction and raw image data body space when promptly being projected as regular projection, directly adopted the image frame data in the raw image data body, generation MIP image;
A02 otherwise, be projected as the free face projection, continue described steps A.
Also comprise the steps a0 before the described steps A 0:
Whether a0, judgement carry out full precision MIP processing, carry out following operation:
A01, if carrying out full precision MIP handles, then continue described steps A 0;
A02 otherwise, be the data volume pixel total amount upper limit with the pixel total amount threshold values T that sets, data volume is made corresponding M IP handles.
Described step a02 comprises the steps:
A021, surpass pixel total amount threshold values T when total pixel of data volume, the data body weight is sampled, with pixel total amount threshold values T is that standard is determined scale down, again according to the scale down of scale down specified data body three-dimensional, each direction is done the space according to this scale down and is equidistantly sampled, and the data after the sampling is made MIP handle;
A022 otherwise, directly data volume is made MIP and handles.
Beneficial effect of the present invention is: in the present invention, principle according to the parallel rays projection, the shape that all section/picture frames project the projecting plane is the same, just has certain skew on the position, if the distance between all contiguous slices/picture frames equates, this side-play amount is equidistant, even if the distance between contiguous slices/picture frame is not equal fully, the relevant conversion of side-play amount is also comparatively simple, the present invention utilizes this characteristic exactly, set up mapping template and mapping table thereof that (the) picture frame forms through projective transformation, for each picture frame thereafter, each value of mapping table is made migration processing, form an interim new mapping table, make MIP according to this table and handle.For the image on projecting plane, only need the image and the magnitude relationship between the existing projected image of relatively newer projection just passable, along with constantly refreshing of projected image, the projected image that obtains at last is exactly final MIP image, do not need to keep intermediate image, reduced demand storage space.Therefore, the present invention has adopted this projection-shift method to reduce operand and has reduced the intermediate storage space requirement.
The present invention utilizes the characteristic on projecting plane in the space analysis geometry, the addressing of complexity calculated be mapped to the projecting plane and carry out, avoided extensive repetitive operation, can obtain same computational accuracy, calculation procedure is few, and computing velocity is fast, and the present invention is with respect to prior art, significantly reduced operation time, saved storage space.Projection equation of the present invention and translation projecting plane all are arithmetics, do not use complicated deformation equation etc., so error is dropped to minimumly in computation process, can reflect the organ with maximal density, the profile of tissue truly.
In the present invention, need not pretreatment time, do not need to generate the raw data field, normal vector field and gradient fields, only setting up the mapping template, needing a small amount of time when searching mapping table, later computing only need add side-play amount and get final product, and can satisfy the requirement that large-scale data field calculates in real time.
In the present invention, adding along with the projected image of each picture frame, by the up-to-date comparative result of the maintenance of dynamic refresh, the MIP image is the data volume of a picture frame on the projecting plane all the time, can not increase, the frame number of data volume is many more, and the advantage of this method is obvious more, the present invention calculates maximum intensity projection MIP image by the method for dynamic refresh, has satisfied doctor/user for the real-time of readding sheet, accuracy requirement.
Description of drawings
Fig. 1 is a maximum intensity projection MIP principle schematic;
Fig. 2 separates the algorithm principle schematic for wrong cutting in the prior art;
Fig. 3 separates algorithm processing procedure synoptic diagram for wrong cutting in the prior art;
Fig. 4 separates algorithm treatment scheme synoptic diagram for wrong cutting in the prior art;
Fig. 5 is raw image data body of the present invention and space coordinates synoptic diagram thereof;
Fig. 6 is a basic controlling flow process synoptic diagram of the present invention;
Fig. 7 is projection coordinate of the present invention conversion synoptic diagram;
Fig. 8 is the embodiment of the invention 1 a concrete control flow synoptic diagram;
Fig. 9 calculates synoptic diagram fast for axial vane surface projection of the present invention;
Figure 10 is the pixel storage mode synoptic diagram of every picture frame of the present invention;
Figure 11 is the embodiment of the invention 2 concrete control flow synoptic diagram.
Embodiment
With embodiment the present invention is described in further detail with reference to the accompanying drawings below:
Embodiment 1:
According to Fig. 5, the raw image data body that constitutes by polylith parallel image frame, picture frame comprise successively the 1st picture frame, the 2nd picture frame ..., and n piece picture frame, form a data volume wide, high, that thickness is respectively W, H, L, adjacent three sides in this data volume space coincide with three-dimensional X-axis, Y-axis and Z axle respectively, and this data volume forms axial vane surface (X-axis * Y-axis), sagittal plane (Y-axis * Z axle) and hat to face (Z axle * X-axis).
As shown in Figure 5, the present invention with the projection of picture frame on the projecting plane wherein as the mapping template, the projection of other picture frame on the projecting plane is benchmark with the mapping template, handle through translation, overlapping with the mapping template, carry out medical image maximal density MIP projection process, because it is the same that all images frame projects the shape on projecting plane, just have certain skew on the position, the distance in the present embodiment between all adjacent image frames equates that this side-play amount is equidistant.Even if the distance between the adjacent image frame is not equal fully, the relevant conversion of side-play amount is also comparatively simple, can determine according to the distance separately between the adjacent image frame.
As shown in Figure 6, basic controlling flow process of the present invention is as follows:
I. as shown in Figure 7, set up the plane parameter of projecting plane P, projecting plane P is vertical with projecting direction, sets up corresponding M IP image two dimension (x axle-y axle) coordinate parameters in the P of projecting plane, with the 1st projection of picture frame on the P of projecting plane as the mapping template.
II. produce mapping table, mapping table reflects the pairing two-dimentional MIP image coordinate location of the three-dimensional space position of data volume in the 1st picture frame.
III. according to mapping table, each point of the 1st picture frame is made medical image maximal density MIP handle.
IV. for follow-up picture frame, be benchmark, each value in the mapping table is made migration processing, form interim new mapping table with the mapping template.
V. according to new mapping table follow-up picture frame is made corresponding medical image maximal density MIP and handle,, produce final MIP image MIP image dynamic refresh.
As shown in Figure 8, concrete control flow of the present invention is as follows:
1. judge the axial anastomose property between projecting direction and the raw image data body, carry out following operation:
1A, if projecting direction and raw image data body space is axially identical, when promptly being projected as regular projection, directly adopt the image frame data in the raw image data body, be divided into following three kinds of situations:
1AA, axial vane surface projection, as shown in Figure 5, projecting direction and Z axle coincide, because the raw image data body is with a frame one frame, store along the mode of Z-direction, therefore basic comparative approach is to get a picture frame with back one frame compares at every turn, carries out the MIP processing.In order to accelerate comparison speed, can adopt the method for the key frame in the video coding to calculate fast, promptly for example shown in Figure 9, when calculating the axial vane surface projection for the first time, getting 16 continuous frames is one group.The pixel value of each relatively present frame and next frame correspondence position writes down current comparative result, and makes comparisons successively with the picture frame of back, finally obtains the maximal value of this 16 frame consecutive image, and the MIP result of these 16 picture frames is formed a key frame.If a total m piece picture frame then can be divided into the k=m/16 group.After the axial vane surface project time, if view field covers a certain group in the k group data, then this group data do not need to compare frame by frame, but get final product with the MIP data participation computing of key frame.Continue following steps 13.
1AB, sagittal plane projection, as shown in Figure 5, projecting direction and X-axis are coincide, in the same picture frame, the basic order that pixel is deposited as shown in figure 10, the location mode of the every frame of image is with row major, so when this direction projection, light is consistent with the direction that each pixel order is arranged, so, only need read in proper order the pixel value of each point and value thereafter once relatively after, what obtain is exactly required MIP value.Take to handle with MIP like the axial vane surface projection-type, just along the X-direction grouping, but not Z-direction in this processing procedure, is made inner linear interpolation processing for the data that generate.Continue following steps 13.
1AC, hat are to the face projection, as shown in Figure 5, projecting direction and Y-axis are coincide, this projecting direction is perpendicular to the X-Z plane, so a picture frame of each structure is a data line of getting a picture frame successively, get the delegation of the correspondence position of back one frame then, and the like to last one deck, one group of data this mode can being formed are thought a picture frame, can take out the picture frame of short transverse like this, these images be done and above similar grouping, in this processing procedure, make inner linear interpolation processing for the data that generate, finish MIP and handle.Continue following steps 13.
1B otherwise, be projected as the free face projection, continue following steps 2.
2. as shown in Figure 7, determine the parameter of projecting plane P, establishing the space projection plane equation is ax+by+cz+d=0, and this projecting plane P is vertical with projecting direction, on the P of projecting plane, choose x direction of principal axis and the y direction of principal axis of two perpendicular directions, make their normal vector be respectively (a as the MIP image
1, b
1, c
1), (a
2, b
2, c
2).
3. get the longest straight line in the raw image data body as the length of side of MIP image, for this regular hexahedron of raw image data body, the longest straight line wherein is a diagonal line, and promptly wide, the height of MIP image are
Wherein W, H, L are respectively wide, the height and the thickness of raw image data body, like this, can make the MIP image possess generality, prevent the image disappearance.
4. get the center of the intersection point of central point on the projecting plane of raw image data body as final MIP image.As shown in Figure 5, raw image data body center point P (x
0, y
0, z
0)=(Width/2, Height/2, Length/2), then, as shown in Figure 7, P (x
0, y
0, z
0) intersection point on the P of projecting plane is P ' (x
0', y
0', z
0'), with P ' (x
0', y
0', z
0') as the center of final MIP image.
5. calculate the 1st three dimensional space coordinate behind the spot projection in the picture frame, according to the intersection point computing formula of putting the plane, the some P (x in the space
0, y
0, z
0) to plane π: the intersection point Q of ax+by+cz+d=0 vertical line (coordinate formula z) is for x, y:
Because the projecting plane just can not influence final result as long as direction is consistent, so can do some reasonable simplification:
If initial point is crossed on the projecting plane, d=0 then.The projecting plane equation simplification is ax+by+c=0.In addition, can be with the coefficient normalization of projecting plane equation, i.e. a
2+ b
2+ c
2=1.
Some P (x in the data volume space
0, y
0, z
0) real space position P ' (x on the P of projecting plane
0', y
0', z
0') calculate by following formula of reduction:
Order:
FTemp＝x
0×a+y
0×b+z
0×c
Then:
x
0′＝x
0-FTemp×a，y
0′＝y
0-FTemp×b，z
0′＝z
0-FTemp×c
P ' (the x that is calculated
0', y
0', z
0') be the some P (x in the space
0, y
0, z
0) three dimensional space coordinate in the P of projecting plane.
6. be the two-dimensional coordinate in the MIP image with the conversion of the volume coordinate after the projection, as shown in Figure 7, calculate P among the P of projecting plane according to following formula
a' (x
a', y
a', z
a') two-dimensional coordinate Position (X in the MIP image
a, Y
a),
Make DeltX=x
a'-x
0', DeltY=y
a'-y
0', DeltZ=z
a'-z
0',
TmpX＝DeltX×a
1+DeltY×b
1+DeltZ×c
1，
TmpY＝DeltX×a
2+DeltY×b
2+DeltZ×c
2，
X then
a=(int) (MIPPlaneSize/2)+TmpX+K
0
Y
a＝(int)(MIPPlaneSize/2)+TmpY+K
0
Wherein, MIPPlaneSize is the width on MIP projecting plane, that is,
At above-mentioned X
a, Y
aCalculating in inserted correction coefficient K
0, to prevent putting in order missing by a mile of computing generation owing to asking, in the present embodiment, K
0Value is 0.5.
Like this, just set up any 1 the pixel Position (X in the MIP image in the 1st picture frame
a, Y
a) between mapping relations one by one.
7. calculate successively in the 1st picture frame the two-dimensional coordinate in the somewhat pairing MIP image, and be stored in the mapping table.
8. according to mapping table, each point of the 1st picture frame is made medical image maximal density MIP handle.
9. for follow-up picture frame, with the mapping template is benchmark, each value in the mapping table is made migration processing, form interim new mapping table,, just certain skew is arranged along Z-direction because other every picture frame is the same with the 1st wide height of picture frame, so can prove from mathematics, since the 2nd picture frame, if certain pixel through the two-dimensional coordinate Position after conversion (X, Y) with the 1st picture frame in the pixel of correspondence position through the two-dimensional coordinate Position (X after the conversion
1, Y
1) differ a Δ x and Δ y, so for all pixels in the same picture frame, this Δ x and Δ y remain unchanged.In other words, the conversion of each picture frame later on being equivalent to that the transformation results of the 1st picture frame has been done individual translation handles.Therefore, we the Δ x and the Δ y that need note in each picture frame get final product.Since the 2nd picture frame, (0,0) locational pixel of each picture frame and (0,0) locational pixel of the 1st picture frame are made comparisons, calculate offset x and Δ y.
If the distance between two width of cloth images is SliceDistan ce, order
tmpSliceDis?tan?ce＝SliceDis?tan?ce×c，
DeltX＝-tmpSliceDis?tan?ce×a，
DeltY＝-tmpSliceDis?tan?ce×b，
DeltZ＝SliceDis?tan?ce-tmpSliceDis?tan?ce×c，
Then:
Δx＝DeltX×a
1+DeltY×b
1+DeltZ×c
1
Δy＝DeltX×a
2+DeltY×b
2+DeltZ×c
2
Since the 2nd picture frame, offset x and the Δ y with respect to the 1st picture frame notes every picture frame.The projection coordinate of other pixel in same picture frame only need add that these two side-play amounts are just passable.Like this, just form interim new mapping table.
10. according to new mapping table follow-up picture frame is made corresponding medical image maximal density MIP and handle, MIP image dynamic refresh.Circular flow step 9 is to step 10, and the 2nd picture frame-n piece picture frame finished medical image maximal density MIP processing in to the raw image data body.
11. the MIP image is carried out the interpolation operation, because between the layer that is written into image and the layer is spaced, and this is at interval common greater than the distance between neighbor in the picture frame, therefore, when carrying out oblique projection, the situation of space or some region projection deficiency might appear.Addressed this problem two kinds of methods: the one, the forward direction interpolation, promptly blank map is as layer data; The 2nd, back to interpolation, i.e. MIP image after the interpolation projection.Better from the later a kind of mode of the angle that improves speed, the present invention adopts gridiron pattern translation interpolating method to solve this problem preferably after obtaining final MIP image.
Gridiron pattern translation interpolating method is actually a kind of interpolation method.From the principle of projection-translation algorithm as can be seen, after the figure layer that inserts in each figure lamellar spacing was made MIP and handled, its position relation was equivalent to that original MIP perspective view has been done a translation and handles.In order to save operation time and space, the present invention does not make linear interpolation to original MIP perspective view, and directly makes the translation overlap-add procedure.Suppose that the figure number of plies that each figure lamellar spacing domestic demand is inserted is M
0, then to be M at least
0Inferior translation overlap-add procedure, the same peaked mode of getting that adopts during stack.
Consider a kind of opposite extreme situations below: when original MIP perspective view is linear, during along interpolation direction (current this direction is the vertical direction of lines) interpolation data, can form the wire space, this is because the straight line in the image is actually discrete form, and discontinuous data, for this reason, need qualification is made in the interpolation path.With a point is example, supposes that its initial position is P (x
0, y
0), need pointwise to be incorporated into P (x
1, y
2) position, we choose the straight line of this point-to-point transmission as the interpolation path, but must walk the gridiron pattern mode, and direct of travel promptly can only be level or vertical, and can not do oblique advancing.Can guarantee that so final MIP image the cavity can not occur.
12. then, adopt the method for 3 * 3 Gauss's smothing filterings further to improve picture quality, like this, the present invention has obtained the MIP image of an optimization by step 11 and step 12, makes the present invention more be tending towards perfect.
13. produce final MIP image.
In fact, three kinds of situations of regular projection among the above-mentioned steps 1A also can be seen the special case of free face projection as, and in other words, three kinds of situations of regular projection also can adopt the disposal route of free face projection.
Embodiment 2:
In the above-described embodiments, embodied full precision MIP treatment scheme to the raw image data body, application for reality, for satisfying real-time Display projector result, the present invention may need to adapt to the change of projecting direction, must will produce sluggishness above certain rise time and feel strict qualification of rise time work of MIP image.
As shown in figure 11, the concrete control flow of another kind of the present invention is as follows:
1) judge that whether carrying out full precision MIP handles, and carries out following operation:
1a) if carrying out full precision MIP handles, then according to full precision MIP flow processing, that is, and according to step 1-step 13 among the embodiment 1, until step 11);
1b) otherwise, be standard according to the pixel total amount threshold values T (in the present embodiment, pixel total amount threshold values T is set at 4M) that sets, carry out following operation:
1b1) the total pixel when data volume surpasses pixel total amount threshold values T, the data body weight is sampled, with pixel total amount threshold values T is that standard is determined scale down, again according to the scale down of scale down specified data body three-dimensional, each direction is done the space according to this scale down and is equidistantly sampled, and continues following steps 2).
1b2) otherwise, continue following steps 2).
2) determine the parameter on projecting plane, choose x direction of principal axis and the y direction of principal axis of two perpendicular directions as the MIP image.
3) length of side of body diagonal length of fetching data as the MIP image.
4) get the center of final MIP image.
5) three dimensional space coordinate behind the spot projection in the 1st picture frame of calculating.
6) conversion of the volume coordinate after the projection is the two-dimensional coordinate in the MIP image.
7) calculate successively in the 1st picture frame the two-dimensional coordinate in the somewhat pairing MIP image, and be stored in the mapping table.
8), each point of the 1st picture frame is made medical image maximal density MIP handle according to mapping table.
9) for follow-up picture frame, be benchmark with the mapping template, each value in the mapping table is made migration processing, form interim new mapping table.
10) according to new mapping table follow-up picture frame is made corresponding medical image maximal density MIP and handle, MIP image dynamic refresh.
11) the MIP image is carried out the interpolation operation.
12) produce final MIP image.
Concrete operation method in the present embodiment and embodiment 1 are described same or similar, and this repeats no more.Present embodiment, is made MIP to the data volume of the pixel total amount threshold values T that resamples/be no more than and is handled as the data volume pixel total amount upper limit by pixel total amount threshold values T, and the time can be controlled in 60 milliseconds, can satisfy the needs of user's true-time operation.Show after tested, the MIP that remakes once full precision when the user stops to change projecting direction handles (aforesaid flow process), and the user is when dragging with mouse, and display speed does not have delay substantially, image detail keeps also more complete, makes the present invention can not produce sluggish sense when showing.
Experimental result proves that the present invention can realize maximum intensity projection MIP fast, and picture quality reaches clinical requirement fully, the operational outfit computing machine adopts the P4/2.66G dominant frequency, and the 512M internal memory is for the image of 300 frame 256*256, under full accuracy modes, reconstruction speed is 281 a milliseconds/frame; Under the non-full accuracy modes of present embodiment, do test with the image of 690 frame 512*512, reconstruction speed is 47 a milliseconds/frame, and this speed can not reduce because of the increase of picture number in theory, can satisfy the requirement of real-time demonstration fully.
In the above-described embodiments, for ease of explanation, the size of all images frame in (original image) data volume is all consistent, in actual applications, even the not of uniform size of picture frame causes, for the image on projecting plane, only need the image and the magnitude relationship between the existing projected image of relatively newer projection just passable, can with all projected images regular be onesize, this can not need to pay creative work to those skilled in the art can implement, and repeats no more herein.
The present invention proposes the image frame data body is converted to the projecting plane two-dimensional coordinate from three dimensional space coordinate conversion method, in the above-described embodiments, for ease of explanation, picture frame adopts the hexahedron of rule, obviously, this conversion method of the present invention is equally applicable to other Processing Algorithm relevant with projection, as many plane projections MPR, curved surface projection CPR or the like, this can not need to pay creative work to those skilled in the art and can implement.
In sum, although ultimate principle of the present invention, method, flow process are specifically set forth by the foregoing description, under the prerequisite that does not break away from main idea of the present invention, according to above-described inspiration, those of ordinary skills can not need to pay creative work can implement multiple conversion/alternative form or combination, repeats no more herein.
Claims (6)
1. maximum density projection generating method for medical image, it is characterized in that: it adopts following steps:
A, for the raw image data body that constitutes by polylith parallel image frame, with the projection of picture frame on the projecting plane wherein as the mapping template; Described steps A comprises the steps:
A1, set up the plane parameter on projecting plane, set up projecting plane corresponding M IP image coordinate parameter, with the 1st projection of picture frame on the projecting plane as the mapping template; Described steps A 1 comprises the steps:
A11, determine the parameter on projecting plane, on the projecting plane, choose x direction of principal axis and the y direction of principal axis of two perpendicular directions as the MIP image;
A12, get the longest straight line in the raw image data body as the length of side of MIP image;
A13, the center of getting the intersection point conduct final MIP image of central point on the projecting plane of raw image data body;
A2, generation mapping table, the pairing two-dimentional MIP image coordinate location of three-dimensional space position of data volume in the described mapping table reflection picture frame, described steps A 2 comprises the steps:
Three dimensional space coordinate behind the spot projection in A21, the calculating chart picture frame;
A22, the volume coordinate after projection conversion is the two-dimensional coordinate in the MIP image;
A23, successively calculate in the 1st picture frame the two-dimensional coordinate in the somewhat pairing MIP image, and be stored in the mapping table;
B, the projection of other picture frame on the projecting plane are benchmark with described mapping template, handle through translation, and be overlapping with described mapping template, carries out the medical image maximum intensity projection and handle, and described step B comprises the steps:
B1, according to described mapping table, each point of the 1st picture frame is made medical image maximal density MIP handles;
B2, for follow-up picture frame, with the mapping template be benchmark, each value in the described mapping table is made migration processing, form interim new mapping table,
B3, the new mapping table of basis are made corresponding medical image maximal density MIP to follow-up picture frame and are handled, to MIP image dynamic refresh, produce final MIP image, among the described step B3, produce before the final MIP image, described MIP image is carried out the interpolation operation, and gridiron pattern translation interpolating method is adopted in described interpolation operation.
2. maximum density projection generating method for medical image according to claim 1 is characterized in that: in the described steps A 22, in the conversion of the two-dimensional coordinate in the MIP image, insert correction coefficient K0.
3. maximum density projection generating method for medical image according to claim 1 is characterized in that: described MIP image is further carried out the operation of Gauss's smothing filtering.
4. according to any described maximum density projection generating method for medical image among the claim 1-2, it is characterized in that: also comprise the steps A0 before the described steps A:
Axial anastomose property between A0, judgement projecting direction and the raw image data body, carry out following operation:
A01, if axially coincideing of projecting direction and raw image data body space when promptly being projected as regular projection, directly adopted the image frame data in the raw image data body, generation MIP image;
A02 otherwise, be projected as the free face projection, continue described steps A.
5. maximum density projection generating method for medical image according to claim 4 is characterized in that: also comprise the steps a0 before the described steps A 0:
Whether a0, judgement carry out full precision MIP processing, carry out following operation:
A01, if carrying out full precision MIP handles, then continue described steps A 0;
A02 otherwise, be the data volume pixel total amount upper limit with the pixel total amount threshold values T that sets, data volume is made corresponding M IP handles.
6. maximum density projection generating method for medical image according to claim 5 is characterized in that: described step a02 comprises the steps:
A021, surpass pixel total amount threshold values T when total pixel of data volume, the data body weight is sampled, with pixel total amount threshold values T is that standard is determined scale down, again according to the scale down of scale down specified data body three-dimensional, each direction is done the space according to this scale down and is equidistantly sampled, and the data after the sampling is made MIP handle;
A022 otherwise, directly data volume is made MIP and handles.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2006100634318A CN101174331B (en) | 2006-11-01 | 2006-11-01 | Maximum density projection generating method for medical image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2006100634318A CN101174331B (en) | 2006-11-01 | 2006-11-01 | Maximum density projection generating method for medical image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN101174331A CN101174331A (en) | 2008-05-07 |
| CN101174331B true CN101174331B (en) | 2011-07-27 |
Family
ID=39422823
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN2006100634318A Active CN101174331B (en) | 2006-11-01 | 2006-11-01 | Maximum density projection generating method for medical image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN101174331B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5675056B2 (en) * | 2008-07-04 | 2015-02-25 | 株式会社東芝 | X-ray imaging apparatus and image processing apparatus |
| CN104796649B (en) * | 2014-01-21 | 2017-12-26 | 北京炬力北方微电子有限公司 | A kind of method and device of tripleplane |
| CN108734705A (en) * | 2018-05-17 | 2018-11-02 | 杭州电子科技大学 | Digital galactophore fault image calcification clusters automatic testing method based on deep learning |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1995019167A1 (en) * | 1994-01-14 | 1995-07-20 | Mallinckrodt Medical, Inc. | Therapeutic treatment for inhibiting vascular restenosis |
| CN1183154A (en) * | 1996-02-06 | 1998-05-27 | 索尼计算机娱乐公司 | Apparatus and method for drawing |
| CN1395904A (en) * | 2001-06-21 | 2003-02-12 | Ge医疗系统环球技术有限公司 | Magnetic resonance inaging system |
| CN1605088A (en) * | 2001-12-20 | 2005-04-06 | 皇家飞利浦电子股份有限公司 | Computer graphics system and method for rendering an image for display |
| CN1688168A (en) * | 2005-06-24 | 2005-10-26 | 清华大学 | Visual difference calculating method for video image in multi-vision point system |
-
2006
- 2006-11-01 CN CN2006100634318A patent/CN101174331B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1995019167A1 (en) * | 1994-01-14 | 1995-07-20 | Mallinckrodt Medical, Inc. | Therapeutic treatment for inhibiting vascular restenosis |
| CN1183154A (en) * | 1996-02-06 | 1998-05-27 | 索尼计算机娱乐公司 | Apparatus and method for drawing |
| CN1395904A (en) * | 2001-06-21 | 2003-02-12 | Ge医疗系统环球技术有限公司 | Magnetic resonance inaging system |
| CN1605088A (en) * | 2001-12-20 | 2005-04-06 | 皇家飞利浦电子股份有限公司 | Computer graphics system and method for rendering an image for display |
| CN1688168A (en) * | 2005-06-24 | 2005-10-26 | 清华大学 | Visual difference calculating method for video image in multi-vision point system |
Non-Patent Citations (2)
| Title |
|---|
| 诸葛斌，冯焕清等.基于有序体数据的最大密度投影算法.中国科学技术大学学报33 4.2003,33(4),全文. |
| 诸葛斌,冯焕清等.基于有序体数据的最大密度投影算法.中国科学技术大学学报33 4.2003,33(4),全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN101174331A (en) | 2008-05-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US4737921A (en) | Three dimensional medical image display system | |
| US10346504B2 (en) | 3D modelling of bodies | |
| CN103106685B (en) | A kind of abdominal organs three-dimensional visualization method based on GPU | |
| CN110060337B (en) | Carotid artery ultrasonic scanning three-dimensional reconstruction method and system | |
| US7586501B2 (en) | Simultaneous projection of multi-branched vessels and their context on a single image | |
| CN103971396B (en) | OpenGL ES (open graphics library for embedded system) implementation method for ray casting algorithm under ARM+GPU (advanced RISC machine+graphic processing unit) heterogeneous architecture | |
| EP1775685B1 (en) | Information processing device and program | |
| CN102525662A (en) | Three-dimensional visual tissue organ operation navigation method and system | |
| US20170161936A1 (en) | Mapping 3d to 2d images | |
| CN102663818A (en) | Method and device for establishing three-dimensional craniomaxillofacial morphology model | |
| Soferman et al. | Advanced graphics behind medical virtual reality: evolution of algorithms, hardware, and software interfaces | |
| CN101290684B (en) | Three-dimensional ultrasound pattern rapid plotting method and apparatus | |
| CN101174331B (en) | Maximum density projection generating method for medical image | |
| CN100476876C (en) | Method for computer-assisted rebuilding heart mitral annulus | |
| CN103699716A (en) | Virtual organ display method driven by personalized three-dimensional medical image | |
| CN101996415B (en) | Three-dimensional modeling method for eyeball | |
| US20070064002A1 (en) | Method and apparatus for providing efficient space leaping using a neighbor guided emptiness map in octree traversal for a fast ray casting algorithm | |
| Guan et al. | Innovative volume rendering using 3D texture mapping | |
| US20090322785A1 (en) | Method, apparatus, system and computer-readable medium for interactive shape manipulation | |
| CN104200511B (en) | Multi-resolution volume rendering method based on intra-block interpolation | |
| EP2266457B1 (en) | Intermediate image generating method, device, and program | |
| Kaufman et al. | A survey of architectures for volume rendering | |
| Subramanyan et al. | A procedure to average 3D anatomical structures | |
| JP5245811B2 (en) | Voxel array visualization device | |
| Yoo et al. | Interactive 3D Medical Visualisation: A parallel approach to surface rendering 3D medical data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| ASS | Succession or assignment of patent right |
Owner name: SHENZHEN LANDWIND NETWORK CO., LTD. Free format text: FORMER OWNER: LANYUN INDUSTRY CO., LTD.;SHENZHEN CITY Effective date: 20120810 |
|
| C41 | Transfer of patent application or patent right or utility model | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20120810 Address after: Futian District King Road Shenzhen city Guangdong province 518034 green garden building room 410 (only office) Patentee after: Shenzhen Landwind Network Co., Ltd. Address before: Futian District King Road Shenzhen city Guangdong province 518034 No. 81 country garden E building 6 floor Patentee before: Landwind Co., Ltd. |
|
| C56 | Change in the name or address of the patentee |
Owner name: LANWANG TECHNOLOGY CO., LTD. Free format text: FORMER NAME: SHENZHEN LANDWIND NETWORK CO., LTD. |
|
| CP03 | Change of name, title or address |
Address after: Nanshan District Xueyuan Road in Shenzhen city of Guangdong province 518000 No. 1001 Nanshan Chi Park A3 building 6 floor Patentee after: LANWON SCIENCE & TECHNOLOGY CO., LTD. Address before: Futian District King Road Shenzhen city Guangdong province 518034 green garden building room 410 (only office) Patentee before: Shenzhen Landwind Network Co., Ltd. |
|
| C56 | Change in the name or address of the patentee |
Owner name: LANWON TECHNOLOGY CO., LTD. Free format text: FORMER NAME: LANWANG TECHNOLOGY CO., LTD. |
|
| CP01 | Change in the name or title of a patent holder |
Address after: Nanshan District Xueyuan Road in Shenzhen city of Guangdong province 518000 No. 1001 Nanshan Chi Park A3 building 6 floor Patentee after: LANWON TECHNOLOGY CO., LTD. Address before: Nanshan District Xueyuan Road in Shenzhen city of Guangdong province 518000 No. 1001 Nanshan Chi Park A3 building 6 floor Patentee before: LANWON SCIENCE & TECHNOLOGY CO., LTD. |