CN102782701B - Patient moving vector is extracted in mark position from radioscopic image - Google Patents

Abstract

A kind of method and system of the patient moving for determining in image.The method includes obtaining image based on the view data produced by scanning device during scanning.Described image comprises at least three mark being assumed to have rigidity or semi-rigid construction.The position of each the actual measurement having in the first dimension and the second dimension on the detector panel of scanning device in described at least three mark.The method farther includes to determine reference three-dimensional position for each in described at least three mark, and defines for being described in relation, the geometric parameter of scanning device and the equation of patient moving between reference three-dimensional position and the actual position measured.The method finally includes: solve these equations in number, think in each the reference three-dimensional position in described at least three mark and described at least three mark corresponding each actual measurement position between the image that takes in of difference, draw the six component motion vectors for describing patient moving.

Description

Patient moving vector is extracted in mark position from radioscopic image

Related application

Present patent application on November 25th, 2009 submit to, entitled " landmark identification in radioscopic image and process (MARKER IDENTIFICATION AND PROCESSING IN X-RAY IMAGES) " Copending U.S. Patent Application No.12/626, the part continuation application of 197, the full content of above application all joins in the application by quoting.This patent Application further relates to that entitled " the patient moving vector using mark position from radioscopic image to be extracted is revised and rebuilds X Ray image (CORRECTING AND RECONSTRUCTING X-RAY IMAGES USING PATIENT MOTION VECTORS EXTRACTED FROM MARKER PO SITIONS IN X-RAY IMAGES) " application No.12/700, 032, entitled " method (the METHOD FOR ACCURATE of the accurate sub-pixel location of the mark on radioscopic image SUB-PIXEL LOCALIZATION OF MARKERS ON X-RAY IMAGES) " application No.12/626,218, entitled " method (the METHOD FOR X-RAY MARKER positioned for X-ray mark in the 3 d space when existing mobile LOCALIZATION IN 3DSPACE IN THE PRESENCE OF MOTION) " application No.12/626,268, Yi Jibiao Entitled " for following the tracks of method (the METHOD FOR TRACKING X-RAY of X-ray mark in CT projection image sequence MARKERS IN SERIAL CT PROJECTION IMAGES) " application No.12/626,310.

Copyright notice

A part of disclosure of patent document comprises material protected by copyright.Patent and trademark office is occurred at it When patent document or record, copyright owner does not oppose that any one in patent documentation or patent disclosure is replicated again Existing, but all retain all copyrights at other any aspects.

Technical field

Embodiments of the invention relate to image procossing.Specifically, embodiments of the invention provide for from radioscopic image In mark position extract patient moving vector, and use the motion for patient's and outside other plan of these vectors Revise the method and system of the image after rebuilding.

Background technology

Modern imaging system obtains many images during single process or scanning.The process being known as rebuilding is used to carry out group Close image, represent, to make, one or more images that the finite width through patient " is cut into slices ".Because for medical treatment or dentistry When process is prepared or in order to suitably calibrate imaging device, often gained image is analyzed, it is essential that Image is the clearest.But, during scanning, patient moves the error or fuzzy that can cause in image.

It is usually present two kinds of methods controlling the impact of patient moving.A kind of method is physically to limit patient.But, thing Reason limits and is often not enough to prevent patient moving, this is because patient can not keep static (such as, young child or due to thing The reason limit) or because operator use patient to limit system inadequately.Even if it addition, carefully carry out patient location and Limit, some patient movings also can occur.

The method of the another kind of impact for controlling patient moving is to measure it when patient moving occurs, and It is modified for it.The existing mechanism moved for detecting patient includes based on laser (optics) and ultrasonic (sound) Range-based searching device (finder) or one group of stereoscope (stereoscopic) camera are with detection patient position change.But, this Detector, range-based searching device or camera used in a little systems adds complexity and the cost of system.Such as, great majority Rotary stand in scanning circumstance makes it difficult to be placed into detector, range-based searching device or camera during stand rotates not Can be on the appropriate location that some moment occurs obstruction.It addition, hardware installation to be added with stand answering of detection process Polygamy.

When Computed tomography (" CT ") imaging, original (raw) view data is reconstructed into final image and relies on specific The fact that physical arrangement changes its location one by one in the most describable mode projection picture.Projection is run through for describing The image of image sequence, structure U and V position has for describing to image (image-to-image) equation of locus formula The parameter of the geometry of stand and the parameter of expression structure physical positioning in three-dimensional (" 3D ") space.If parameter Fixed, then equation describes track exactly.But, if parameter changes (such as, owing to patient transports in scanning process Dynamic), then actual path is different from expected trajectory.Mathematical relationship between measured position and the desired location of structure of structure It it is the basis extracting the motion vector representing patient moving.

After detection motion, there is various ways with for the Motion correction image detected.A kind of mode is mobile X Source detector is to offset the impact of motion.This is main method for image stabilization in camera.Another kind of method is The information about the motion detected is used after obtaining view data, it to be modified.This can include that deletion is laid equal stress on The new movement obtained owing to detecting and destroyed image.But, this method cannot process continuous print patient during scanning Motion (that is, destroying all parts or the motion of the overwhelming majority of image), and add obtain required for image time The area of a room.It addition, the reacquisition of view data is infeasible in some form of single rotation sweep program.

The another kind of method being modified for the motion detected in the picture is to attempt to compensate for detected motion And the mode making image be revealed as seeming generation of not moving translates, scales and rotate image.This method can effectively be located in Reason patient'sTranslationMotion (translational motion) but it can not process patient's completelyRotateMotion.Such as, suffer from Two points on person can be different, but after patient rotates, these 2 can be revealed as a point on gained image.Figure As translation or scaling can not be modified for this situation, at least it is generally not capable of.

Another revises the method for the motion detected as described in U.S. Patent application No.2004/0136590 Make like that to rebuild grid (grid) warpage (warp).Make to rebuild grid warpage and can process deformable motion, such as, exist The motion occurred during heart movement.Therefore, this method uses the information about the motion detected to revise and rebuilds grid, with Make itself it will be assumed that deformation account for.Although this method can process the motion of broad range, but form grid and incite somebody to action It is complicated that view data matches grid.Therefore, this method can be time-consuming, and by various types of artifacts (artifact) impact.

Summary of the invention

Embodiments of the invention provide for determining six component patient position (moving) vectors for each projection picture System and method.A kind of method includes, occurs without patient moving, then obtain three on each projection picture or more Each U and V position in multiple X-ray benchmark (fiducial) mark, and determine to explain the actual location of mark and mark The patient position skew required for change between the expection location of will.Based on described change, the method output describes projection The translation of patient and six component vector of rotary motion in Xiang.The method allows to exist with about 0.2 millimeter or preferably degree of accuracy Patient moving is followed the tracks of continuously during whole scanning.

Embodiments of the invention also provide for for being attached in process of reconstruction about the information of the patient moving detected Method and system, its use the information to compensate motion.A kind of method includes obtaining six components about each projection picture Each vector is also combined, to generate about respectively by position error vector (that is, patient moving vector) with static bench alignment information The projection matrix of individual projection picture.Described projection matrix describes each X-ray is how to go through three dimensional object from x-ray source Enter two-dimensional x-ray detector.Described process of reconstruction uses projection picture and described projection matrix to rebuild image and to this Image is modified for the patient moving detected.

In a particular embodiment, it is provided that a kind of method of patient moving for determining in image.The method includes, Image is obtained by view data produced by scanning device based on during scanning.Described image include being assumed having rigidity or At least three mark of semi-rigid construction, wherein, each in described at least three mark has in the first dimension and the second dimension The position of the actual measurement on the detector panel of scanning device on degree.Then by electronic processing unit determine about described at least Each reference three-dimensional position in three marks.Electronic processing unit then define or be provided for describe reference three-dimensional position and Relation, the geometric parameter of scanning device and the equation of patient moving between the actual position measured.Finally, electron process These equations are solved by unit in number, thus show that for described image six component motion describing patient moving are vowed Amount, described six component motion vectors are to each the reference three-dimensional position in described at least three mark and described at least three mark In will corresponding each actual measurement position between difference take in.

In another specific embodiment, based on the patient moving correction detected and rebuild multiple projection picture.The method Including, at computer, obtain the multiple projection pictures produced by scanning device.Each in the plurality of projection picture include to Few three marks, and each in described at least three mark have in the first dimension and the second dimension at scanning device Measured position on detector panel and measured three-dimensional position.Each based in the plurality of projection picture In at least three mark, by electronic processing unit determine about each the site error in the plurality of projection picture to Amount.Patient moving in described position error vector definition projection picture.By about each in the plurality of projection picture Each position error vector with and scanning device be associated geometric parameter combination, with draw for the plurality of projection as in Each projection matrix.The plurality of projection picture and each the projection matrix that is used in the plurality of projection picture are carried Supplying or be supplied to image reconstruction process, this process generates the image for the revised reconstruction of patient moving.

Other aspects of the present invention will be by taking in detailed description and accompanying drawing and become obvious.

Accompanying drawing explanation

Fig. 1 shows the cone-beam computed tomography system according to an embodiment.

Fig. 2 diagrammatically illustrates parts and the geometrical relationship thereof of the cone-beam computed tomography system of Fig. 1.

Fig. 3 diagrammatically illustrates the computer of the Fig. 1 according to an embodiment.

Fig. 4 shows the image comprising mark.

Fig. 5 shows the flow chart that the image after rebuilding according to an embodiment provides the process of Motion correction.

Fig. 6 graphically illustrates the three class rotary motions of the patient of the horizontal plane relative to detector panel and suffers from The two class translational motions of person.

Fig. 7 graphically illustrates the three class rotary motions of the patient of the vertical plane relative to detector panel and suffers from The two class translational motions of person.

Fig. 8 shows three parts (three-part) mistake for determining patient moving vector according to an embodiment The flow chart of journey.

Fig. 9 shows the flow chart of the image correcting method according to an embodiment.

Detailed description of the invention

Before any embodiments of the invention are explained in detail, it is to be understood that the present invention is not limited to when it is applied That illustrated in explained below or structure illustrated in the accompanying drawings details and the layout of parts.The present invention can have Other embodiments, and can be practiced or carried out in every way.

Also should give it is noted that may utilize multiple devices based on hardware and software and the portion of multiple different structure Part realizes the present invention.Additionally, as described in chapter subsequently, concrete configuration shown in figure is intended to illustrate the reality of the present invention Execute example.Alternative arrangements is also possible.

Modern X-ray process uses one or more detectors and non-film.Detector electronically detects and quantifies to arrive Reach the amount (that is, the density of X-ray) of the X-ray of detector.Using detector, setting up X ray computer tomography art (" CT ") is System, it makes x-ray source rotate around patient, and utilizes the single wide strip in the side relative with x-ray source of patient (single-wide strip) detector element, electronically detection gained X-ray.X-ray source, detector and permission The combination of the frame for movement that x-ray source rotates is known as " CT stand ".Collect Autonomous test for all different X-ray positions The data of device, are then combined it during being known as rebuilding.Graphical representation after combination is single through patient's Finite width " is cut into slices ", and the image density at the most each point represents the X-ray density of the tissue at specific physical positioning.Rebuild Process uses the fixed geometry of CT stand and the changing angle of the x-ray source-detector panel combination relative to patient (that is, θ) processes collected projection picture.

By repeating said process while patient mobile between image acquisition or scanning, x-ray source or both, Multiple section can be obtained.Such as, the platform of mobile supporting patient and the most also movable gantry, then replace data slicer and Produce " helical (helix) ".It addition, (such as, detector strip or the size of ring or width are increased to multirow from single-line detector Reach 256 row), it is allowed to more data is more quickly gathered.Additionally, replace detector strip with bigger two-dimensional detector, Then obtain the whole detector panel image at each x-ray source position, rather than only wall scroll data.To quantitatively reaching 600 or the collection of more these image be known as projection picture.Each projection graphical representation is along with x-ray source and detection Device synchronously rotate around patient and from different perspectives or angle the X-ray snapshot to patient.Because cone type X-ray bundle needs Cover two-dimensional detector, so this CT imaging is known as cone beam (" CB ") CT imaging.Fig. 1 diagrammatically illustrates according to this The CB CT imaging system 10 of a bright embodiment, and Fig. 2 shows parts and the geometrical relationship thereof of CB CT imaging system 10 And parameter.

CB CT imaging system 10 includes scanning device 12 and computer 14.Scanning device 12 includes CT stand 13, and it includes that X penetrates Line source 16 and detector 18.X-ray source 16 and detector 18 are rotating on microscope carrier 20 at the alignment of opposite each other, this rotation microscope carrier 20 X-ray source 16 and detector 18 is made to move around patient 22.Patient 22 is supported by seat 24.Imaging system shown in Fig. 1 10 is dental imaging system.Therefore, it is placed in bracket (rest) 26 during patient 22 is sitting in seat 24 and by his or her lower jaw. When stand 13 has rotated the scanning of the head of patient, bracket 26 keeps the head geo-stationary of patient.

View data is exported to computer 14 by scanning device 12.Pictorial data representation is examined by detector 18 during scanning The intensity level of the X-ray measured.Computer 14 is connected to control station 30, and it includes display 32 and one or more input And/or output device, such as keyboard 34.User uses control station 30 with mutual with computer 14.Such as, user can use Control station 30 is to ask image or other data from computer 14.The information asked is provided to control station 30 by computer 14, And this information is shown on the display 32 by control station 30, print this information to printer (not shown), and/or by this letter Breath preserves to computer-readable memory module (not shown).

Fig. 3 diagrammatically illustrates the computer 14 of Fig. 1 according to an embodiment of the invention.Computer 14 includes defeated Enter/output interface 40, electronic processing unit (" EPU ") 42 and one or more memory module, all drive such as by disk The computer readable disk of dynamic device (not shown) access, random access memory (" RAM ") module 44, read only memory (" ROM ") module 46 or a combination thereof.Input/output interface 40 receives view data from scanning device 12, and this view data is carried Supply EPU 42.

In certain embodiments, complete stand rotates time-consuming about 8 ~ 40 seconds.During this period, patient may move or CT stand may unexpectedly move, thus causes the fuzzy of gained image.Typical image resolution ratio is the amount of 0.25 millimeter Level.Therefore, the patient moving of this same order frequently results in image blurring, and large-scale patient moves and may make gained figure As becoming unacceptable for the purpose of its desired outpatient service.Even if it addition, when obscuring excessive, obscuring and also can Picture quality is caused to be generally reduced.Such as, stand vibration can cause picture quality that is fuzzy and that reduce.

Computer 14 is moved for patient by the movement of rigidity or semi-rigid object in tracing figure picture and carries out image Revise.Such as, not having under the ideal conditions that patient moves wherein, the rigid objects being imaged rotates around patient along with stand And change the location in the two dimension of projection picture in clearly defined mode.Patient moves (or the stand outside expection moves) and leads Cause the deviation between the expection location of objects in images.Therefore, by measuring the inclined of clearly defined object and its desired locations From, can measure and revise patient's (and stand outside expection) amount of movement.Specifically, as described below, if image sequence Row exist at least three object, then can combine these measured objects and expect that with it deviation positioned is to determine that patient transports Dynamic vector, it can apply to image and is modified to move for patient.

In order to ensure there being the clearly defined rigid objects of desired quantity in the picture, can be before the scan by base Fiducial mark will (such as, three or more) is placed on patient.These marks generally comprise lead or steel BB, and it is comparatively dense and can prevent Or restriction X-ray transparent.But, these marks can also be made up of other materials, and is configured to generation during scanning In projection picture with relatively high ratio visible other form or shape.Each mark or multiple mark can be located at binding agent Between Ceng, and this binding agent can coat to patient, to guarantee that mark does not moves during program.Although it is it addition, the most multiple Miscellaneous but it also may the anatomic landmark thing (landmark) within Shi Yonging is as mark, rather than the mark of applications.

Mark is placed on patient so that each visual field or the image that are generated by CB CT imaging system all include at least Three marks.Such as, patient can place seven to nine marks, with guarantee at least three mark in each image surface, with Reduce position measurement noise, and allow the statistical combination of result.In certain embodiments, on patient, mark is uniformly spaced Open, to have the spatial distribution of maximum about patient, and avoid interference image interpretation (interpretation).

After mark is placed on patient, patient is scanned, the head of such as patient, and would indicate that projection picture The gained view data of sequence is transmitted to computer 14.Fig. 4 shows the exemplary projected image 100 including eight marks 102. Should give it is noted that " mark 102 " be the expression of mark rather than the mark of reality (that is, mark 102 is not actual bb or class As device).Use term mark to refer to actual tag and the expression of mark or image is used with those skilled in the art Universal grammar syntax consistent.Therefore, in following written description,NotAlways to actual mark and the table of mark Make between showing and clearly distinguishing.

EPU 42 receives view data and processes information by performing one or more application softwaries or module.Apply soft Part or module are stored in memory module, such as ROM module 46 or computer readable disk or can be by EPU 42 by driving The medium moving device or external interface and access.As it is shown on figure 3, memory module stores mark processing module 50, motion vector carries Delivery block 52, projection matrix computing module 54 and reconstruction module 56.Fig. 5 shows that method 60, the method 60 show through this The data stream of four modules.As it is shown in figure 5, EPU 42 obtains the projection picture from input/output interface 40, and perform mark Processing module 50, with identify index point on each projection picture, identify each identified mark vertically and horizontally position, Each position and dimension thereof are distributed to suitable physical token and uses image location data to estimate the benchmark of each mark Three-dimensional (" 3D ") physical positioning (step 62 and 64).It is No.026212-9015-00, title in the application attorney docket distributed For " landmark identification in radioscopic image and process (MARKER IDENTIFICATION AND PROCESSING IN X-RAY IMAGES) copending application " describes the various embodiments of this landmark identification and processing procedure.

As it is shown in figure 5, obtain after the marker point list of each projection picture at EPU 42, EPU 42 performs patient Motion vector draws module 52, to draw six component patient moving vectors (step 66 and the 67) (ginseng for each projection picture See Fig. 6-8).EPU 42 then performs projection matrix computing module 54, with several based on six component patient moving vectors and stand What parameter and determine the projection matrix (step 68) after the Motion correction of each projection picture (participation Fig. 9).Then EPU 42 perform to rebuild module 54, to generate Motion correction volumetric image (step 69) based on the projection matrixes after Motion correction (seeing Fig. 9).The each several part of detailed hereafter method 60.

Patient moving vector draws

As it is shown in figure 5, the first step of method 60 includes the marker position information obtaining each mark in projection picture (step 62).It is " the mark knowledge in radioscopic image No.026212-9015-00, entitled in the application attorney docket distributed Not and process (MARKER IDENTIFICATION AND PROCESSING IN X-RAY IMAGES) " in copending application Disclose the various embodiments of marker position information for determining each mark in projection picture.Described marker position information U(row including each mark on each projection picture) and V(row) location.U and the V location of each mark provides and scanning device 2D position (that is, the U of the relevant mark of detector panel_M1、V_M1、U_M2、V_M2、U_M3And V_M3).These measured values are to there is trouble The position of the actual measurement of the mark detected when person moves.As shown in Figures 6 and 7, the first dimension U is circular with the rotation of CT stand Tangent.Second dimension V is the plane that the rotation axis with CT stand is paralleled.

In step 64, method 60 determines the benchmark 3D physical set indicated for each based on U and the V positional information indicated Position (that is, X, Y and Z).Each motion indicated is judged by the reference position indicated.Reference position is arbitrary, and can be each The mode of kind is defined.In one embodiment, it is defined as the position of mark when scanning starts.This position can be with various sides Formula draws, such as finds X, Y and the Z coordinate being used for generating the theory locus of the track of the mark observed by best fit.When making During with three marks, for these origin reference locations indicated by being labeled as X₁、Y₁、Z₁、X₂、Y₂、Z₂、X₃、Y₃And Z₃Nine coordinates It is worth (that is, three coordinates that three marks are multiplied by each) to be defined.

Following equation mathematical way describes location X_i、Y_iAnd Z_iThe benchmark physics 3D point at place is how to project to spy Determine stand rotational value θ_jIn the 2D detector plane at place:

$U_{\mathrm{ij}}=\frac{\mathrm{DSD}*(Y_i\cos {\mathrm{\θ}}_j-X_i\sin {\mathrm{\θ}}_j)}{(\mathrm{DSO}+Y_i\sin {\mathrm{\θ}}_j+X_i\cos {\mathrm{\θ}}_j)}$ (equation 1)

$V_{\mathrm{ij}}=\frac{\mathrm{DSD}*Z_i}{(\mathrm{DSO}+Y_i\sin {\mathrm{\θ}}_j+X_i\cos {\mathrm{\θ}}_j)}$ (equation 2)

Subscript i represents that tag number, subscript j represent projection numbering.DSO be from x-ray source to stand center of rotation away from From, and DSD is to the distance of detector panel from x-ray source.Variable θ_jIt it is current gantry rotational angle.Such as, when stand source When being positioned at the head dead astern of object, θ_jIt is zero.Fig. 2 generally illustrates parameter DSO, DSD and θ.

Then in step 66, method 60 determine containing for specific projection image six kinematic variableses (that is, Δ X, Δ Y, Δ Z, α, β and γ) patient moving vector.Generally, when the object scanned is considered rigid body, this step determines description, and this is swept Retouch the six component motion vectors for projection picture of the overall movement of object (such as, the head of patient).Each motion vector All consider the observed position (U and V) of each mark in projection picture and the difference of the prediction location of each mark corresponding, institute State prediction location suppose preferable stage motion and each mark known, fixing (such as, non-moving) physics X, Y and Z origin reference location.

Such as, in this step, computer receives the list object for each mark as input, described object column tabular Go out the positional information for the mark in each projection picture.If there is at least three mark in each projection picture, Then computer is based on the logo object list of mark included in this image, and exports the expression for each projection picture The six component motion vectors for the motion detected by this image.Three expression patients in six components tie up at X, Y and Z Translation displacement (shift) on degree or motion (that is, variable Δ X, Δ Y and Δ Z), the other three representation in components patient is about X, Y Rotation displacement or motion (that is, variable α, β and γ) with Z axis.

As it is shown in figure 5, do not rely on other projection pictures, determine patient moving vector (step for each projection picture 67).Therefore, step 66 and 67 end product be that one group of motion of a large amount of vectors comprising the quantity equal to projection picture is vowed Amount.Such as, 600 projection pictures are studied generation is had the motion vector set of 600 Independent Vector.But, if one Individual or multiple projection pictures do not comprise three or more marks, then the actual quantity of produced vector can be less, and this can hamper Method 60 is hindered to draw movable information.

Should give it is noted that in certain embodiments, motion vector is the estimation of the accurate motion vector of projection picture Value.Such as, during scanning, each mark experiences different change at its X, Y with Z location.Method 60 uses each projection In Xiang, the different motion of each mark is to generate six component motion vectors.If there is three marks, then can use six sides Formula (that is, each U and the V equation in three marks) and six unknown quantitys (that is, the component of motion vector) are come really Fixed answer accurately or motion vector.But, if image exists the mark of more than three, then, compared with unknown quantity, exist More equations, and this problem is overdetermination (overspecified).It addition, mark measurement can not be always accurate due to noise Really measure.Accordingly, because problem is overdetermination and there is measurement noise, so being to use by motion vector determined by method 60 Estimated value in the accurate motion vector of specific projection image.

Return to step 66, various equations and formula can be defined, to determine patient moving vector.Such as, as There is not patient moving in fruit, then for the X of special sign_i、Y_iAnd Z_iBeing all constant for all projection pictures, this implies U and V Merely due to change θ and change.But, when there is patient moving, X_i、Y_iAnd Z_iChange due to the positioning effects of motion.Cause This, when there is patient moving, each coordinate is more suitably represented by following equation:

X_i=X_iref+x_ij Y_i=Y_iref+y_ij Z_i=Z_iref+z_ij

Wherein, x_ij、y_ijAnd z_ijRepresent special sign i position with it at specific projection j being associated with patient moving The difference of origin reference location.Plug these values in equation 1 and 2 illustrated above, obtain following revised equation:

$U_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*\left(\right(Y_{\mathrm{refi}}+y_{\mathrm{ij}})*\cos {\mathrm{\θ}}_j-(X_{\mathrm{refi}}+x_{\mathrm{ij}})*\sin {\mathrm{\θ}}_j)}{(\mathrm{DSO}+(X_{\mathrm{refi}}+x_{\mathrm{ij}})*\cos {\mathrm{\θ}}_j+(Y_{\mathrm{refi}}+y_{\mathrm{ij}})*\sin {\mathrm{\θ}}_{\text{j}})}$ (equation 3)

$V_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*(Z_{\mathrm{refi}}+z_{\mathrm{ij}})}{(\mathrm{DSO}+(X_{\mathrm{refi}}+x_{\mathrm{ij}})*\cos {\mathrm{\θ}}_j+(Y_{\mathrm{refi}}+y_{\mathrm{ij}})*\sin {\mathrm{\θ}}_j)}$ (equation 4)

During when using three marks, at specific θ_jPlace, exist six equations (for three mark in each Individual U and V) and nine unknown quantity (each x in three marks_ij、y_ijAnd z_ij), it is not having the feelings of additional information Under condition, can not solve.

But, if increasing in all marks of supposition are all fixed on rigid body or this restriction on rigid body, the most any specific mark The motion of will can be by its origin reference location and by six transformation parameters (that is, motion vector components) Δ X, Δ Y, Δ Z, α, β The motion of the whole object provided with γ completely describes.Variable Δ X, Δ Y and Δ Z define the translation of patient in projection picture Move, and variable α, β and γ represent the rotary motion of patient in projection picture, as shown in Figures 6 and 7.Therefore, by increasing rigid body Limiting, so far only exist six unknown quantitys and six equations, this makes problem to solve.

For further reduced equation and formula, the infrastructural frame for object can fixing (i.e., relatively from benchmark In ground) framework changes into the rack frame of the rotation of benchmark.Rotating frame conversion can be would be attached to convert with patient moving Combination, to provide overall movement to convert.Lower column matrix provides this combination:

|Γ_transformed|=|A_overall||Γ_ref|=|θ_rotation||A_{patient_motion}||Γ_ref| (equation 5)

Wherein, | Γ | is the position vector of the object comprising component X, Y and Z, | А_{patient_motion}| it is to comprise and Δ X, Δ The motion transform matrices of the component that Y, Δ Z, α, β and γ translation and rotation are associated, | θ_rotation| it is for specific projection The spin matrix being associated with stand rotation of picture.The place of equation of the equation 3 and 4 that this hint is above will is that

$U_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*Y_{\mathrm{transformes\; j}}}{\mathrm{DSO}+X_{\mathrm{transformed}}}$ (equation 6)

$V_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*Z_{\mathrm{transformed\; j}}}{\mathrm{DSO}+X_{\mathrm{transformed}}}$ (equation 7)

Wherein:

X_transformedj=f₁(X_refi,Y_refi,Z_refi, Δ X, Δ Y, Δ Z, α, beta, gamma, θ) (equation 8)

Y_transformedj=f₂(X_refi,Y_refi,Z_refi, Δ X, Δ Y, Δ Z, α, beta, gamma, θ) (equation 9)

Z_transformedj=f₃(X_refi,Y_refi,Z_refi, Δ X, Δ Y, Δ Z, α, beta, gamma, θ) (equation 10)

And f₁、f₂And f₃Defined by equation of transformation (that is, equation 5) above.

Therefore, " fixing " parameter DSD, DSO, X are given_ref1、Y_ref1、Z_ref1、X_ref2、Y_ref2、Z_ref2、X_ref3、Y_ref3And Z_ref3And the U at specific θ_MEAS1、V_MEAS1、U_MEAS2、V_MEAS2、U_MEAS3, and V_MEAS3Measured value, method 60 user's formula 5,6 and 7 draw at this θ, and kinematic parameter Δ X at each other θ, Δ Y, Δ Z, α, β in the projection image sequence And γ.

It should be understand that, shown equation 6 and 7 itself is nonlinear.By 5,8,9 and 10 tables of equation The motion transform shown has sine (sine) and cosine (cosine) component that the rotating part with conversion is associated.And, by In each middle intrinsic division, thus solve expression nonlinear operation while equation 6 and 7.All of these factors taken together is the darkest Show and cannot draw Exact Solutions, and iterative will be the most time-consuming.But, by using some to simplify, it may be determined that can The solution of row.Fig. 8 shows the flow process of three Part Methods 70 for determining patient moving according to an embodiment of the invention Figure.When EPU 42 perform patient moving vector draw module time, the EPU 42 of computer 14 perform method 70.Method 70 uses Two simplification, are all based on actual patient moving premise relatively small compared with stand rotation and carry out.

First simplifies by the SIN function (sine) during rotation transformation and cosine function (cosine) are applied Taylor Series expansion, and only consider each in Section 1 carry out lienarized equation formula.Effectively, it means that cos (x) replaces with 1 And sin (x) x replaces.This supposition is effective for the rotational component of little value.But, as described below, method 70 is permissible Use and allow this to limit the alternative manner relaxed.

Second simplification be advantageously employed V equation to the change of α, β and Δ Z the most sensitive and to Δ X, Δ Y and γ Change less sensitive this point.On the contrary, U equation to the change of Δ X, Δ Y and γ the most sensitive and to α, β and the change of Δ Z Change less sensitive.It means that for little quantity of motion, general problem can carry out in addition letter by PROBLEM DECOMPOSITION becomes two parts Change.Specifically, it is possible, firstly, to the V equation after linearisation while α, β and Δ Z is solved, then can be by The value for α, β and Δ Z that previously used V equation draws, as constant, is come linearisation while Δ X, Δ Y and γ After U equation solve.

As it has been described above, use iteration can improve the degree of accuracy solved.Iteration (step 76) comprises overlappingly to U equation (step 72) and V equation (step 74) solve.In the solution of each equation, existing about iteration is how to increase degree of accuracy Two components.First, V each solve all use by U previous solve and draw for Δ X, Δ Y and γ nearest , the value that degree of accuracy gradually steps up.Equally, U each solve all use by V previous solve and draw about α, β and Δ Z Nearest value (step 74).Second component is the use of residual error (residual).In order to accomplish this point, by U and V equation Higher-order (nonlinear component) be considered as pseudo-constant, to allow linear formula.Owing to suffering from for normal stage motion Person's motion may be considered little, it is possible to makes this supposition.In U or V equation, any one each solves (step 71,72 or 74), after, pseudo-constant is updated to the value that its degree of accuracy gradually steps up.First for V equation solves (step 71), before performing to solve, there is no information to set for pseudo-constant or Δ X, the value of Δ Y and γ.Performing this step Before, these parameters are set to " zero ".

Provide below and solve V equation and the further detail below of U equation.

Solve V equation

In order to solve V equation, it is noted that component α is defined as the rotation around X-axis, β is defined as the rotation around Y-axis, and γ defines For rotation about the z axis.Accordingly, the matrix form of the rotation represented with α is as follows:

$\mathrm{rot\α}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\mathrm{\α}\right)& -\sin \left(\mathrm{\α}\right)\\ 0& \sin \left(\mathrm{\α}\right)& \cos \left(\mathrm{\α}\right)\end{array}\right]$

To this matrix application Taylor expansion, make the SIN function in this matrix and cosine function linearisation.

$\cos \left(\mathrm{\α}\right)=1-\frac{{\mathrm{\α}}^2}{2!}+\frac{{\mathrm{\α}}^4}{4!}+...$

$\sin \left(\mathrm{\α}\right)=\mathrm{\α}-\frac{{\mathrm{\α}}^3}{3!}+\frac{{\mathrm{\α}}^5}{5!}+...$

For low-angle, Section 11 and α represents the reasonable approximation of cos (α) and sin (α) respectively.Alternatively, can pass through Degree of accuracy is improved by the estimated value correction Section 1 of higher item.But, in the case of there is no existing understanding, these higher-orders The best estimate of item is zero.But, once the linear approximation being used for sin (α) and cos (α) is used for drawing the approximation of α, Then this α value may be used for being modified to the approximation of higher order term.In other words, can be by being given containing the expression formula of following two Cos (α) and each the expression of sin (α): linear (in α) item and non-linear residual error, this non-linear residual error is in any iteration phase Between can be seen as being constant.

Such as, providing below about cos (α) and the mathematical notation of sin (α), wherein cos (α) replaces with c1, and sin (α) replace with s1 with simpler and clearer.

cos(α)=c1=1-h1

h1=1-c1

sin(α)=s1=a1-g1

g1=a1-s1

Wherein, h1 represents that cos (α) residual error and g1 represent the residual error for sin (α).Suffix numeral 1 is used for indicating rotation alpha.

Rotation β, suffix numeral 3 expression rotate γ to use suffix numeral 2 to represent, and suffix numeral 4 represents that stand rotates θ, Lower column matrix defines each rotation paid close attention to.

$\mathrm{rot\α}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1-h1& -a1+g1\\ 0& a1-g1& 1-h1\end{array}\right]$ (equation 11)

$\mathrm{rot\β}=\left[\begin{array}{ccc}1-h2& 0& -a2+g2\\ 0& 1& 0\\ a2-g2& 0& 1-h2\end{array}\right]$ (equation 12)

$\mathrm{rot\γ}=\left[\begin{array}{ccc}c3& -s3& 0\\ s3& c3& 0\\ 0& 0& 1\end{array}\right]$ (equation 13)

$\mathrm{rot\θ}=\left[\begin{array}{ccc}c4& -s4& 0\\ s4& c4& 0\\ 0& 0& 1\end{array}\right]$ (equation 14)

Note, a1=α and a2=β, and owing to it is to need two variablees solved to be decomposed.On the contrary, γ and θ Solve period at V equation and be accordingly to be regarded as preset parameter due to each in these variablees thus be not decomposed.

For the translational component of motion, 3 component reference positions of following vector representation special sign.

${\mathrm{position}}_{\mathrm{original}}=\left[\begin{array}{c}X0\\ Y0\\ Z0\end{array}\right]$

The mark conversion from a position to another position can relate to up to 3 displacement (translation) and 3 rotations Transfer position.The accurate rank of these displacements and translation is arbitrary, but must be defined for concrete formula.Such as, motion Conversion can be defined as first comprising three-component translation displacement, is followed by α rotation, β rotation and γ and rotates, and is finally θ (platform Frame) rotate.This definition is in following derivation, but equivalent is also suitable for by other definition.

Can be by only increasing the moving displacement (displacement) (XD, YD and ZD) the corresponding benchmark position to mark Put, define new X, the Y for special sign and Z location.But, in this V equation solves, X and Y location are assumed to often Amount, it is possible to by using two symbol XS and YS to simplify.Use these symbols, after the translation of following vector representation mark Position (be applicable to translate mobile and non-rotating after position).

${\mathrm{position}}_{\mathrm{tramslated}}=\left[\begin{array}{c}X0+\mathrm{XD}\\ Y0+\mathrm{YD}\\ Z0+\mathrm{ZD}\end{array}\right]=\left[\begin{array}{c}\mathrm{XS}\\ \mathrm{YS}\\ Z0+\mathrm{ZD}\end{array}\right]$

Location after being applied with rotation, after so far mark is in its complete moving displacement.This new position is claimed For the position vector after conversion, and be given in the following equation.

$\left[\begin{array}{c}{X}_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\end{array}\right]=\mathrm{rot\θ}\×\mathrm{rot\γ}\×\mathrm{rot\β}\×\cot \mathrm{\α}\×\left[\begin{array}{c}\mathrm{XS}\\ \mathrm{YS}\\ Z0+\mathrm{ZD}\end{array}\right]$ (equation 15)

Now, V equation can solve.Such as, it is presented above V equation as equation 6, repeats conduct side here Formula 16:

$V=\frac{\mathrm{DSD}*Z_{\mathrm{transformed}}}{\mathrm{DSO}-X_{\mathrm{transformed}}}$ (equation 16)

Equation 16 represents between the mark location after the V of a special sign and the motion transform of this special sign Relation.If combined with equation 11,12,13,14 and 15 by equation 16, then gained equation comprises three solved afterwards Individual variable: a1, a2 and ZD(correspond to α, β and Δ Z).

0=-DSD*YS*g1-DSD*YS*h2*a1+DSD*h2*h1*Z0+V1*c4*c3*g2*Z0+V1 * XS*c4*c3- V1*YS*c4*s3+DSD*h2*h1*ZD-V1*YS*s4*c3+V1*s4*s3*a2*ZD-V1*DSO-DSD*h1*Z0-DSD*h2* ZD+V1*YS*s4*s3*a2*a1-V1*XS*s4*s3-DSD*XS*g2+DSD*YS*a1-DSD*h2*Z0+DSD*XS*a2-DSD* h1*ZD+DSD*YS*h2*g1+V1*c4*c3*g2*ZD+V1*YS*a2*c4*c3*g1-V1*s4*s3*g2*ZD-V1*s4*c3* g1*ZD+V1*YS*c4*c3*g2*a1-V1*YS*a2*s4*s3*g1-V1*YS*s4*s3*g2*a1-V1*YS*c4*c3*g2*g1 +DSD*Z0+DSD*ZD-V1*c4*s3*g1*ZD-V1*s4*c3*g1*Z0+V1*YS*s4*s3*g2*g1-V1*s4*s3*g2* Z0-V1*c4*s3*g1*Z0+V1*c4*s3*a1*Z0+V1*c4*s3*a1*ZD+V1*s4*c3*a1*Z0+V1*s4*c3*a1* ZD-V1*c4*c3*g2*h1*Z0-V1*s4*s3*a2*h1*Z0-V1*s4*s3*a2*h1*ZD+V1*XS*s4*s3*h2-V1* XS*c4*c3*h2+V1*YS*c4*s3*h1-V1*c4*c3*a2*Z0-V1*c4*c3*a2*ZD+V1*s4*s3*a2*Z0-V1* YS*c4*c3*a2*a1-V1*c4*c3*g2*h1*ZD+V1*s4*s3*g2*h1*Z0+V1*s4*s3*g2*h1*ZD+V1*YS* s4*c3*h1+V1*c4*c3*a2*h1*Z0+V1*c4*c3*a2*h1*ZD

One such equation is all existed for each mark.

In order to solve aforesaid equation, it is necessarily arranged to following form:

A1*A+a2*B+ZD*C+D=0 (equation 17)

This will be by will comprise only one of three variablees a1, a2 or ZD and not under not higher than first power (power) state These subexpressions containing any SIN function or cosine function residual error g1, g2, h1 or h2 are divided into one group to complete.These three groups sons Expression formula becomes coefficient A, B and C.Every other item is gathered the most together in " constant " item D.

Once carrying out this classification, (each mark has one just to obtain the m system of linear equations of the form shown by equation 17 Individual equation).If there is 3 marks, then draw solution by simple matrix inversion.If more than the mark of three Will, then can use singular value decomposition (" SVD ") to draw solution.SVD provides the side that order is associated with each equation effectively Value that difference minimizes, for a1, a2 and ZD.

For first time iteration, constitute the higher-order combination comprising a1, a2 and/or ZD or containing SIN function or cosine The subexpression of the D of function residual error g1, g2 or h1 or h2 is assumed to zero.But, in subsequent iterations, it is provided that for this The estimated value of the improvement of a little higher-order items.

Solve U equation

For representing Δ X, variable X D of Δ Y and γ, the solving of YD and a3 respectively, it is possible to use provide and provide specific change Change the equation of the U position of rear X and Y location, and to realize in the way of solving V equation all fours:

$U=\frac{\mathrm{DSD}*Y_{\mathrm{transformed}}}{\mathrm{DSO}-X_{\mathrm{transformed}}}$

In the case, equation 4 is the most identical, and equation 11,12,13 and 15 becomes:

$\mathrm{rot\α}=\left[\begin{array}{ccc}1& 0& 0\\ 0& c1& -s1\\ 0& s1& c1\end{array}\right]$ (equation 18)

$\mathrm{rot\β}=\left[\begin{array}{ccc}c2& 0& -s2\\ 0& 1& 0\\ s2& 0& c2\end{array}\right]$ (equation 19)

$\mathrm{rot\γ}=\left[\begin{array}{ccc}1-h3& -s3+g3& 0\\ s3-g3& 1-h3& 0\\ 0& 0& 1\end{array}\right]$ (equation 20)

$\left[\begin{array}{c}{X}_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\end{array}\right]=\mathrm{rot\θ}\×\mathrm{rot\γ}\×\mathrm{rot\β}\×\cot \mathrm{\α}\×\left[\begin{array}{c}X0+\mathrm{XD}\\ Y0+\mathrm{YD}\\ \mathrm{ZS}\end{array}\right]$ (equation 21)

After replacing it, it is necessary to store gained equation, to provide the equation of following form:

$\mathrm{XS}*E+\mathrm{YD}*F+a3*G+H=0$ (equation 22)

Furthermore, one of these equations will be there are for each mark.Matrix inversion can be utilized or utilize singular value to divide Solution solves these equations.

Although equation 17 and 22 can be solved by any order, but it it is preferred embodiment first solving equation formula Solving equation formula 22 after 17, then by the two equation iteration two or more times.Have been found that twice iteration can carry For the most useful estimated value for motion, but preferably there is iteration more times.Only five times iteration are utilized to be obtained with essence Exactness (more preferable than six decimal digits).Because when N is typically less than 10, each iteration all comprises and only inverts 3 × N number of square Battle array, it is possible to be performed quickly and totally solve.

Appendix A provides the example of the most above-mentioned method for solving U and V equation Embodiment, it uses Matlab programming language.In Intel(Intel) on Core 2U7600CPU, it is right to be used for performing The calculating time of 24 marks at each place in 300 projections and the most compiled Matlab code of 5 iteration is about Being 0.2 second, this represents that for the approximate solution time of each projection be 7 milliseconds.In this calculating time, major part is divided with singular value Dematrix solves and is associated.

Returning to Fig. 5, after determine six component motion vectors for projection image sequence, EPU 42 examines at step 68 and 69 Rope also performs projection matrix computing module 54 and rebuilds module 56, to generate that revise and after rebuilding projection picture.These steps Rapid detailed description (seeing Fig. 9) in the next section.

It should be noted that in the above-mentioned description to method 70, it is assumed that stand (that is, revolves except expectation in perfect mode Beyond Zhuaning, there is not stand and rotate) rotate around patient.But, can relax this condition, and can amending method 70 to consider plan Outside stage motion.Also can amending method 70, to be applicable to more complicated stand geometry and rotation.

Additionally, in method 70, mark without having fixed relationship each other when being positioned on patient.Between mark Relation can by mark track behavior stagewise (piecewise) analyze and draw, as long as patient moving is the most serious.So And, in certain embodiments, once it being placed, mark is accomplished by mutually maintaining fixed relationship.This means method 70 preferably It is applicable to rigid body, the head of such as people.But, if using more mark, then method 70 can also be used for the less body of rigidity.Example As, as long as three (or more) mark has fixed relationship for specific region, then the region more than can be carried out Evaluate.This design is useful for evaluating the motion of lower jaw and the upper jaw in head discretely when jaw moves.It addition, If three marks are fixed on rigid frame, then can also evaluate the motion of this framework.The method for anatomically can Amoeboid movement, it is useful for such as breathing.

It addition, the mark of method 70 is placed the most critically important for measurement error, this is because mark places impact side The degree of accuracy of the result of method 70.Such as, getting and more open between mark, measurement error will be the least.In one embodiment, around trouble 180 degree of the head of person places mark, and this provides about 0.2 millimeter or more preferable motion vector components detection degree of accuracy.

Because method 70 is each projection picture produces independent motion estimated values, provide additional flexible in this way Property.First, if the mark less than three being detected in specific projection image, then the method can be by adjacent projection picture Motion in this projection picture is carried out interpolation (interpolate).Secondly, the method can use to be had with each projection picture Adjacent projection result is filtered by the data closed, and further provides for the correcting action of this process.

Additionally, because method 70 use imaging data to extract movable information, so the motion detected includes that patient moves Stand outside dynamic and unmodified or plan moves both.Stand outside can using about plan moves Information dynamically revises stage motion.Such as, using method 70, correction can patient scan (such as, if patient by with Instrument and equipment) period carry out, or as use multiple mark anthropometric dummies (phantom) patient scan between calibration A part for journey.

Image correction and reconstruction

Returning to Fig. 5, in projection picture after detection motion, computer 14 is rebuild based on the motion detected and revises this A little images (step 68 and 69).Fig. 9 shows the method 80 performing these steps according to an embodiment of the invention.Work as EPU When 42 execution projection matrix computing modules 54 and reconstruction module 56, the EPU 42 of computer 14 perform image correcting method 80.

Generally, method 80 is by being attached to the information about the motion detected in process of reconstruction for Motion correction Projection picture.Method 80 processes non-variable shape motion detected in the picture, and it is by six component vector (that is, patient moving Vector) specified for each projection picture.As it has been described above, six component vector include three for X, Y and Z translational motion Parameter (that is, Δ X, Δ Y and Δ Z) and for around three parameters (that is, α, β and γ) of X, Y and the rotary motion of Z axis.Each to Amount represent patient moving (such as, in the case of externally measured motion) or patient relative to stand motion (such as, if Movable information deduces from image).Position error vector can represent patient moving correlated error, dynamic stage position error, Or both.

Return to Fig. 9, as the initial step of method 80, obtain CB CT projection picture (step 82).Then, throw for each Shadow image determines position error vector (that is, patient moving vector) (step 84).Can be by including that external detector and image divide All kinds of devices of parser, determine position error vector.In one embodiment, method 80 uses and describes above with reference to Fig. 5-8 Method 60 complete this step.

Then, method 80 is by position error vector and the static calibration parameter or the information (example that describe scanning device geometry As, invariable in scanning process) parameter combination, to draw the projection matrix (step 86) for each projection picture. Parameter for scanning device geometry includes the distance (DSO) from x-ray source to center of rotation, from x-ray source to detector Distance (DSD) and the θ value of each projection picture.Projection matrix describes each X-ray is how to wear from x-ray source traveling Cross patient and specified point (spot) on reaching detector panel.Projection matrix is equal to for operating at 3D computer graphical In on the computer screen virtual 3D volume is converted into the matrix of 2D image.Exemplary projection matrix is as follows:

ξ=A00*X+A01*Y+A02*Z+A03*1.0

ψ=A10*X+A11*Y+A12*Z+A13*1.0

ζ=A20*X+A21*Y+A22*Z+A23*1.0

U=ξ/ζ

V=ψ/ζ

Coordinate X, Y and Z represent the position in the 3 d space of the point on patient, U and V represents that identical point is at 2D detector panel On position, and A00 to A23 represents the component of the projection matrix for giving projection picture.

After generating projection matrix, the projection matrix of projection picture and correspondence is supplied to process of reconstruction by method 80 (88).Described process of reconstruction uses these projection pictures and projection matrix to determine at referred to as rear-projection (back-projection) Period how by each projection as on pixel-map to projection volume.Alternatively, process of reconstruction can use the weight of iteration Build and revise and rebuild image, its use projection matrix with by mapping pixel data to volume data.Process of reconstruction output is rebuild After image, these images are revised based on the motion (data 90) detected.

Even if when not measuring patient moving, method 80 also provides for the mechanism being modified for dynamic stage errors, As long as stage errors is repeatably between scans.Method 80 can also be by using Graphics Processing Unit (" GPU ") Or other special parallel-processing hardware and further speed up.It addition, image correcting method 80 may be adapted to have incomplete fortune The situation of dynamic data.Such as, if the position of mark only starts in scanning and is known at end, then middle projection picture can By utilizing the known position error vector of other projection pictures and to the position error vector being associated with this projection picture Carry out linear interpolation to be revised.

Although patient moving vector being shown that method 60 and image correcting method 80 are retouched already in connection with CB CT imaging State, but these methods can be used for CT, MRI, ultrasonic, the medical imaging of other forms and non-medical imaging, such as, photograph. Such as, image correcting method 80 may be used for the imaging technique wherein evaluated by 3D volume as 2D image sequence, including 3D is ultrasonic and 3D optical homogeneity x-ray tomography art (coherence tomography).

Additionally, although ROM module 46 is shown as storage by Fig. 3 single module (such as, 50,52,54,56 etc.), but Can also combine or separate the module being stored in ROM module 46 or other memorizeies.Such as, in certain embodiments, identical Module perform patient moving vector draw method 60 and image correcting method 80.It addition, ROM module 46 could be included for holding Other modules of row function in addition to the above function.Additionally, in certain embodiments, module can be stored in readable On disk, and it is sent to RAM module 44 and/or ROM module 46 to perform.

Therefore, in addition to other factors, the patient moving that the present invention is provided to determine in image and based on being examined The patient moving correction surveyed the method and system rebuilding image.The various feature and advantage of the present invention will be in appended claims Book illustrates.

Appendix A

1%Load U and V data (n markers x m projections) into

2load′data\test.mat′；

3

4%Load Reference 3D X, Y, and Z positions of each marker (3xn)

5load data\RefLoc.mat；

6

7%Scanner Geometry Parameters

8DSD=716.59；

9DSO=488.111；

10SCANANGLE=360.03；

11

12%Global Parameters

13numMarkers=size (u_mm, 2)；

14numProj=size(u_mm,1)；

15theta=

(-30:-(SCANANGLE/numProj) :-30+ (SCANANGLE/numProj) * numProj+1)) '/ 180.0*pi；

16

17%Initialize motion parameters to zero (first guess)

18CU=zeros (size (u_mm, 1), 3)；

19CV=zeros (size (u_mm, 1), 3)；

20

21%The following starts the main iterative loop

22for iter=1:5

23%*****Solve for V associated parameters (Z, Alpha, and Beta)

24for j=1:numProj

25for i=1:numMarkers

26%The following are " pseudo-constants " in the lineari zed equations

27V1=v_mm (j, i)；

28XS=RefLoc (1, i)+CU (j, 1)；

29YS=RefLoc (2, i)+CU (j, 2)；

30Z0=RefLoc (3, i)；

31ZD=CV (j, 1)；

32a1=CV (j, 2)；

33a2=CV (j, 3)；

34h1=1-cos (CV (j, 2))；

35h2=1-cos (CV (j, 3))；

36g1=CV (j, 2)-sin (CV (j, 2))；

37g2=CV (j, 3)-sin (CV (j, 3))；

38s3=sin (CU (j, 3))；

39c3=cos (CU (j, 3))；

40s4=sin (theta (j))；

41c4=cos (theta (j))；

42

43Alpha_vect (j, i)=DSD*YS+V1* (-s4*c3*Z0+c4*s3*Z0)；

44Beta_vect (j, i)=DSD*XS+V1* (-s4*s3*Z0-c4*c3*Z0)；

45delZ_vect (j, i)=DSD；

46errorV_A (j, i)=(V1*DSO-DSD*Z0+V1* (XS* (-s4*s3-c4*c3)+YS* (c4* S3-s4*c3))) ...

47...

48+ (V1* (XS* (c4*c3*h2+s4*s3*h2) ...

49+YS* (c4*c3*g2*g1-a2*c4*c3*g1-s4*s3*g2*a1-c4*c3*g2*a1+ S4*c3*h1...

50-a2*s4*s3*g1+s4*s3*a2*a1+c4*c3*a2*a1+s4*s3*g2*g1-c4* S3*h1) ...

51+ZD* (c4*c3*a2+s4*c3*a1+s4*s3*g2*h1-s4*c3*g1+c4*c3*g2* H1-c4*c3*g2...

52-s4*s3*g2+c4*s3*g1-s4*s3*a2*h1-c4*s3*a1+s4*s3*a2-c4* C3*a2*h1) ...

53+Z0* (-s4*s3*a2*h1-c4*c3*a2*h1-s4*c3*g1-s4*s3*g2...

54+c4*c3*g2*h1-c4*c3*g2+c4*s3*g1+s4*s3*g2*h1)) ...

55+DSD* (XS*g2+YS* (h2*a1-h2*g1+g1) ...

56+ZD* (h1-h2*h1+h2)+Z 0* (h1-h2*h1+h2)))；

57end

58

59BV=errorV_A (j :) '；

60AV=[delZ_vect (j :)；Alpha_vect (j :)；Beta_vect (j :)] '；

61CV (j :)=(AV BV) '；

62end

63

64%*****Solve for U associated parameters (X, Y, and gamma)

65forj=1:numProj

66for i=1:numMarkers

67%The following are " pseudo-constants " in the linearized equations

68U1=u_mm (j, i)；

69X0=RefLoc (1, i)；

70Y0=RefLoc (2, i)；

71XD=CU (j, 2)；

72YD=CU (j, 2)；

73ZS=RefLoc (3, i)+CV (j, 1)；

74c1=cos (CV (j, 2))；

75s1=sin (CV (j, 2))；

76c2=cos (CV (j, 3))；

77s2=sin (CV (j, 3))；

78a3=CU (j, 3)；

79g3=a3-sin (a3)；

80h3=1-cos (a3)；

81s4=sin (theta (j))；

82c4=cos (theta (j))；

83

84delX_vect (j, i)=(U1*c2*c4-DSD*c2*s4)；

85delY_vect (j, i)=(U1*c1*s4-U1*s1*s2*c4+DSD*s1*s2*s4+DSD*c1* c4)；

86gamma_vect (j, i)=(U1* (X0*c2*s4-Y0*s1*s2*s4-Y0*c1*c4+ZS*s1* C4-ZS*c1*s2*s4) ...

87+DSD* (X0*c2*c4+Y0*c1*s4-Y0*s1*s2*c4-ZS*s1*s4-ZS*c1*s2* c4))；

88errorU_vect (j, i)=U1*DSO-U1* (X0*c2*c4-Y0*s1*s2*c4+Y0*c1*s4- ZS*s1*s4...

89-ZS*c1*s2*c4) ...

90-DSD* (-X0*c2*s4+Y0*c1*c4+Y0*s1*s2*s4+ZS*c1*s2*s4-ZS* S1*c4) ...

91...

92+U1* (XD* (-c2*s4*a3+c2*s4*g3+c2*c4*h3) ...

93+X0* (c2*c4*h3+c2*s4*g3) ...

94+YD* (c1*c4*a3-s2*s1*c4*h3+s2*s1*s4*a3+s2*s1*s4*g3-

C1*c4*g3...

95+c1*s 4*h3) ...

96+Y0* (-s2*s1*s4*g3-s2*s1*c4*h3-c1*c4*g3+c1*s 4*h3) ...

97+ZS* (-s1*s4*h3-s2*c1*c4*h3+s1*c4*g3-s2*c1*s4*g3)) ...

98...

99+DSD* (XD* (c2*c4*g3-c2*c4*a3-c2*s4*h3)+X0* (c2*c4*g3-c2* S4*h3) ...

100+YD* (-c1*s4*a3+c1*c4*h3+c1*s4*g3+s2*s1*c4*a3+s2*s1* S4*h3...

101-s2*s1*c4*g3) ...

102+Y0* (c1*c4*h3+s2*s1*s4*h3-s2*s1*c4*g3+c1*s4*g3) ...

103+ZS* (-s1*s4*g3-s2*c1*c4*g3+s2*c1*s4*h3-s1*c4*h3)) ...

104；

105end

106

107AU=[delX_vect (j :)；DelY_vect (j :)；Gamma_vect (j :)] '；

108BU=errorU_vect (j :) '；

109CU (j :)=(AU BU) '；

110end

111end

112

113figure；Plot ([CV (:, 2), CV (:, 3), CU (:, 3)] * 180/pi)；

114title (' Motion Angles:Alpha, Beta, and Gamma ')；

115

116figure；Plot ([CU (:, 1), CU (:, 2), CV (:, 1)])；

117title (' Motion Translations:X, Y, and Z ')；

Claims (13)

1., for the method determining the patient moving in image, described method is by including with electronic processing unit and storage The imaging system of the computer of the memory module of the motion vector extraction module that can be performed by described electronic processing unit is held OK, described method includes:

At described computer, obtaining image based on the view data produced by scanning device during scanning, described image comprises Being assumed at least three mark with rigidity or semirigid structure, each in described at least three mark has first In dimension and the second dimension, the position of actual measurement on the detector panel of scanning device；

Utilize described electronic processing unit, determine each reference three-dimensional position in described at least three mark；

Utilizing described electronic processing unit, definition is for describing each described reference three-dimensional position in described at least three mark And relation, the geometric parameter of described scanning device and the equation of patient moving between the described actual position measured；And

Utilize described electronic processing unit, described equation is solved, think that described image draws for describing patient's fortune Six dynamic component motion vectors, described six component motion vectors are to each described reference three-dimensional position in described at least three mark The difference put between the position that in described at least three mark, each described reality is measured accordingly takes in.

Method the most according to claim 1, wherein, in described at least three mark, each has reference three-dimensional position X_i、Y_i And Z_iWith the measurement position U on the detector panel of described scanning device in the first dimension and the second dimension_mAnd V_m, wherein, institute State six component motion vectors and include Δ X, Δ Y, Δ Z, α, β and γ, and wherein, described equation is carried out numerical solution to obtain Go out six component motion vectors to include:

A () supposes Δ X, the value of Δ Y and γ,

(b) use Δ X, the assumed value of Δ Y and γ and define the described equation of position in described second dimension wherein it One, calculate the value for Δ Z, α and β, and

(c) use the value of calculation of Δ Z, α and β and define position in described first dimension described equation one of them, Calculate for Δ X, the value of Δ Y and γ,

Δ X, Δ Y, Δ Z represent patient's translation displacement on X, Y, Z-dimension respectively；α, β, γ represent that patient is about X, Y respectively Rotation displacement with Z axis.

Method the most according to claim 2, wherein, carries out numerical solution to show that six component motion are vowed to described equation Amount farther includes: repeat step (b) to (c) at least one times.

Method the most according to claim 2, wherein, carries out numerical solution to show that six component motion are vowed to described equation Amount includes: use the SIN function in the Taylor expansion approximation rotation transformation being made up of the residual error of linear Section 1 and estimation And cosine function (cosine) (sine).

Method the most according to claim 4, wherein, carries out numerical solution to show that six component motion are vowed to described equation Amount farther includes: arrives (c) period repeating step (b) every time, updates the residual error of approximated SIN function and cosine function.

Method the most according to claim 4, wherein, uses described Taylor expansion to approximate the sinusoidal letter in described rotation transformation Number and cosine function include the residual error using the estimation as constant.

Method the most according to claim 1, wherein, the described image obtained at described computer includes projection picture.

8., for determining a system for the patient moving in image, described system includes:

Scanning device；With

Computer, described computer includes electronic processing unit, memory module and input/output interface, wherein, described in deposit Memory modules is configured to the motion vector extraction module that storage is performed by described electronic processing unit；

Wherein, described scanning device is configured to be obtained image based on during scanning by view data produced by described scanning device And sending the image to described computer, described image comprises and is assumed have rigidity or at least the three of semirigid structure Individual mark, each in described at least three mark has in the first dimension and the second dimension, in the inspection of described scanning device The position of the actual measurement on survey device panel,

Wherein, described electronic processing unit is configured to determine that each reference three-dimensional position in described at least three mark Put,

Wherein, described electronic processing unit is configured to definition for describing each described benchmark in described at least three mark Relation, the geometric parameter of described scanning device and the side of patient moving between three-dimensional position and the described actual position measured Formula,

Wherein, described electronic processing unit is configured to solve described equation, thinks that described image draws for retouching Stating six component motion vectors of patient moving, described six component motion vectors are to each described base in described at least three mark In quasi-three-dimensional position and described at least three mark, the difference between the position of each described actual measurement accordingly takes in.

System the most according to claim 8, wherein, each in described at least three mark has reference three-dimensional position X_i、 Y_iAnd Z_iWith the measurement position U on the detector panel of described scanning device in the first dimension and the second dimension_mAnd V_m, its In, described six component motion vectors include that Δ X, Δ Y, Δ Z, α, β and γ, and described electronic processing unit are configured to Following steps carry out numerical solution to draw six component motion vectors to described equation:

A () supposes Δ X, the value of Δ Y and γ,

(b) use Δ X, the assumed value of Δ Y and γ and define the described equation of position in described second dimension wherein it One, calculate the value for Δ Z, α and β, and

(c) use the value of calculation of Δ Z, α and β and define position in described first dimension described equation one of them, Calculate for Δ X, the value of Δ Y and γ,

Δ X, Δ Y, Δ Z represent patient's translation displacement on X, Y, Z-dimension respectively；α, β, γ represent that patient is about X, Y respectively Rotation displacement with Z axis.

System the most according to claim 9, wherein, described electronic processing unit is configured to by repeating step Suddenly (b) to (c) carrys out to carry out described equation numerical solution at least one times to draw six component motion vectors.

11. system according to claim 9, wherein, described electronic processing unit is configured to use by linear one SIN function in the Taylor expansion approximation rotation transformation of the residual error composition of secondary item and estimation and cosine function, come described equation Formula carries out numerical solution to draw six component motion vectors.

12. system according to claim 11, wherein, described electronic processing unit is further configured to by weight every time Multiple step (b) updates approximated SIN function and the residual error of cosine function to (c) period, and described equation is carried out numeral Solve to draw six component motion vectors.

13. systems according to claim 11, wherein, described electronic processing unit is configured to use as constant The residual error of estimation use described Taylor expansion to approximate the SIN function in described rotation transformation and cosine function.