System and Method for Determining Internal Radioactivity and Absorbed Dose Estimates

1. Field of the Invention

The present invention relates to a system and method for determining internal radioactivity and absorbed dose estimates of a patient. More specifically, the present invention determines internal radioactivity uptake in a patient who is given a radiopharmaceutical or other radioactive carrier by using a CT scan assisted matrix inversion (CAMI) method, and the present invention also estimates radiation dose received by a patient by using a Monte Carlo assisted voxel source kernel (MAVSK) method.

2. Background and Related Art

Radiology consists of several major divisions including diagnostic radiology, therapeutic radiology, and special radiologic imaging systems. Whereas the main function of the diagnostic radiologist is to produce and interpret images of internal organs of the body using radiation, therapeutic radiology concerns itself primarily with the treatment of disease by the destruction of diseased tissue with high-energy radiation. Special imaging radiology includes nuclear techniques, computerized tomography (CT), and ultrasound.

In nuclear medicine, radiation is usually internally distributed in the body and arises from the decay of radioactive materials which a patient receives (e.g., by injection, inhalation, or ingestion). An important application is organ imaging, in which the patient receives a gamma-emitting radioisotope that accumulates in an organ (or conversely, the exclusion of an isotope from a particular organ).

Pathological lesions, variations in blood perfusion or anatomical variations will produce an abnormal accumulation (or exclusion) of the isotope in particular organs. Accordingly, mapping the distribution of the isotope by detecting the emitted gamma rays yields an image which provides clues of the presence of pathological lesions.

Mapping the time distribution of the radioisotope image usually requires following the activity over one or several selected areas. Devices for this include the simplest types of radiation detector systems (probes) or gamma cameras which generally use scintillation crystals because they are the most efficient for gamma detection.

The basic configuration of such a probe generally consists of a collimator to limit the field of view, followed by a scintillator, which in turn is followed by a photomultiplier with associated electronics for recording. Output curves are handled in many ways, such as by use of oscilloscope display, scalers, count-rate meters, or strip-chart recorders.

In computerized tomography, a very narrow x-ray beam is passed through the body in a plane perpendicular to the axis of the body. Total attenuation of the beam is then measured by a detector. This measurement is repeated thousands of times with different rays across the body in the same plane to obtain thousands of attenuation measurements across a narrow slice of the body in a plane perpendicular to the axis.

The attenuation of a narrow beam along a path can be expressed as a product of the attenuation through each segment of that path. Since the attenuation is exponential, the total attenuation argument .mu.x can be expressed as the sum of the individual arguments for each segment:

With this concept, one can envision a CT slice as being composed of a number of small volume elements, called voxels. Each measurement of the attenuation of the narrow beam across a slice is thus a measurement of the attenuation of those voxels in the path of that measurement. Accordingly, it is possible to set up a set of linear equations for the solution to the attenuation in each voxel element.

In clinical and experimental radioimmunoimaging and radioimmunotherapy (RIT), knowledge of the absorbed radiation dose in normal tissue and in tumors is very important for an understanding of normal tissue toxicity and tumor dose response. For example, the dose due to monoclonal antibodies (Mabs) is calculated rather than measured from radioactive content of the tissues. Absolute radioactivity distributions are usually calculated from some non-invasive imaging procedures. The gamma camera is a non-invasive method currently available to assess normal organ and tumor biodistribution of radiolabeled antibodies.

Conventional methods of calculating internal radioactivity have various limitations. For example, although some investigators believe that quantitative single photon emission computer tomography (SPECT) will become the method of choice, it is currently difficult to implement and error-prone.

In the opposing view first developed by Thomas et al. (in "In Vivo Quantitation of Lesion Radioactivity Using External Counting Methods," Med. Phys., Volume 3, pp. 253-255 (1976) and "Quantitative External Counting Techniques Enabling Improved Diagnostic and Therapeutic Decisions in Patients With Well-Differentiated Thyroid Cancer," Radiology, Volume 122, pp. 731-737 (1977)), regions of interest (ROIs), or the areas which may be medically significant, are drawn around each organ on both anterior and posterior planar gamma-camera scans. The total radioactivity in the organ is estimated by use of various information including the geometric mean of the anterior and posterior views of the pertinent area, an attenuation correction, a calibration factor of the camera, the radioactive counts per minute from the anterior and posterior views, and a buildup factor B.

More specifically, for a point source in the anterior projection:

From the posterior projection:

where

.mu.=attenuation coefficient (cm.sup.-1);

d=source depth measured from anterior surface (cm);

T=total patient thickness at the source position (cm);

B(x)=buildup factor (a depth (x) dependent parameter to account for phantom scatter. This factor was not used in Thomas' original work.)

C.sub.A =counts per minute (cpm) in anterior gamma camera view;

C.sub.P =cpm in posterior gamma camera view;

A.sub.O =source radioactivity (.mu.Ci);

CF=camera calibration factor (cpm/.mu.Ci), measured with a source having a known amount of radioactivity.

Assuming the depth dependence of B(x) is very small, i.e., B(x).congruent.1, the source activity can then be calculated by taking the geometric mean of C.sub.A and C.sub.P such that: ##EQU1## For a volume source, the idea can be easily extended: ##EQU2## The factor f=2sinh(.mu.l/2)/(.mu.l), where l is the thickness of the source, is due to the volume distribution of the radioactivity. When background activity exists, an ROI can be drawn near the source organ and background counts can be subtracted from the organ ROI counts on a pixel by pixel basis. Notice that only a single source organ is considered in this analysis and the correction for the contribution from background activity is only a first-order approximation.

However, with high background radioactivity and/or overlapping organs, it is difficult to draw ROI's around the organ, and the first-order background correction may not be accurate enough.

When organs are overlapping, it is impossible to separate the activity from different organs. For example, the liver overlaps the right kidney, and the kidney's activity densities are often too low to show up to allow accurate kidney ROI drawing.

Additionally, there is the problem of an organ or tumor being visible in only one view; e.g., the spleen or kidney from the posterior projection.

Furthermore, some organs are not visible in either anterior or posterior views; e.g., the kidneys at later time points after injection of Mabs. In such cases, the geometric mean is zero or close to zero and the opposing view method cannot be applied.

Thus, when the boundaries of organs do not show up very well on gamma camera images, supplemental scans may be necessary. For example, when the liver is imaged following the injection of .sup.99m Tc-sulfur colloid, the lung is imaged with .sup.99m Tc-macroaggregated albumin, or the kidney is imaged with .sup.99m Tc-DTPA. These extra scans are used only for the outline of the normal organs and give extra radiation doses to the patient.

As noted above, radioactivity in the blood pool generally gives relatively high background tissue counts, especially at early time points.

In order to obtain absolute uptake in these cases, the present invention determines absolute uptake using uptake information from the gamma camera images and anatomic data from CT scans.

The absolute uptake may then be used to determine the absorbed dose received by the patient. In many radionuclide procedures, absorbed dose in various organs from sources within the body has to be estimated.

It is essential to determine the radiation dose in radioimmunotherapy where relatively large amounts of radioactivity are labeled onto anti-tumor antibodies and injected into the patient. The radionuclide is then transported to the tumor site(s) by the tumor specific antibody to deliver radiation dose to the tumor. Since a large amount of radioactivity is involved in this procedure, tumor dose as well as normal organ dose can no longer be just a crude estimation. Both absorbed doses have to be calculated in a more precise and systematic manner.

Conventional methods of absorbed dose estimation suffer from various shortcomings. The Medical Internal Radiation Dose (MIRD) Committee of the Society of Nuclear Medicine provides a general framework for the absorbed dose estimation of administered radiopharmaceuticals, and its most recent practical software implementation is the MIRDOSE3 program.

In MIRDOSE3, a uniform distribution of unit cumulated activities within each source organ is assumed. The average absorbed dose to a given target organ from a given source organ is then estimated by dividing the total energy deposited by the target organ mass (S factors).

This S matrix is calculated for several standard anatomic models (e.g., Standard Adult, Standard Adult Female, 15-year-old Phantom and several other models). The product of cumulated activity in a source organ and corresponding S factor gives the dose contribution from source to the target organ. Total target organ dose is then the summation of all source organ absorbed dose contributions.

There are several limitations with MIRDOSE3. First, S factors are calculated based on standard mathematical human phantoms. Actual patient anatomy in terms of organ shape and size as well as inter-organ position might be very different from any available phantom. Thus the S factors for any individual may vary considerably from the standard phantom.

Second, MIRDOSE3 assumes a uniform distribution of activity in the source organs. It has no way to incorporate non-uniform distribution of the source activity when the information becomes available.

Third, output from MIRDOSE gives only an average absorbed dose in the target organs. No spatially varying radiation dose or statistical information of the dose distribution within an organ (e.g., a volume-dose histogram) is available except for red marrow.

Fourth, the dose estimates for beta particles are calculated assuming complete local absorption. No consideration is given to the range of the beta particles; i.e., cross-fire between adjacent organs.

When the radionuclide is used for diagnostic purposes, normal tissue doses are well below tissue radiation tolerance. Therefore, MIRDOSE3 has sufficient accuracy for absorbed dose estimates in diagnostic imaging studies.

One conventional alternative method uses a "real-time" Monte Carlo calculation for radioimmunotherapy dosimetry. The advantage of a Monte Carlo approach is its flexibility and accuracy. Patient-specific information in terms of anatomy and tissue inhomogeneity can all be incorporated into a Monte Carlo simulation. However, a drawback of Monte Carlo is its extensive computational time. Monte Carlo calculations are by nature time-consuming, because a statistically acceptable answer requires numerous pseudo-particles to be generated and their histories followed.

Several methods to do patient-specific dose estimates have been proposed by different groups. One method proposes a "real-time" Monte Carlo calculation for radioimmunotherapy dosimetry (T. K. Johnson and R. L. Vessella, "On the Possibility of `Real-Time` Monte Carlo Calculations for the Estimation of Absorbed Dose in Radioimmunotherapy," Compute. Meth. Progr. Biomed., Vol. 29, pp. 205-210 (1989) and T. K. Johnson and S. B. Colby, "Photon Contribution to Tumor Dose from Considerations of 1311 Radiolabeled Antibody Uptake in Liver, Spleen and Whole Body," Med. Phys., Vol. 20(6), pp. 1667-1674 (1993)).

The MABDOSE code (of T. K. Johnson et al.) currently executes in a "reasonable" amount of time on a CRAY Y-MP/832 supercomputer. For a general institution without large computational power, the Monte Carlo approach is probably not feasible currently. A second problem with a Monte Carlo simulation is that the surface of the organs described in Monte Carlo code is usually required to be in analytical form (e.g., sphere, cylinder, plane or a volume composed by those surfaces). In a clinical situation, patient's organs have an arbitrary geometric shape. Thus it is difficult to directly use the CT 3-D anatomy in any Monte Carlo calculation without using crude approximations.

Another conventional alternative approach uses a point source function (PSF). The PSF is conceptually simple and easy to implement yet it can give patient-specific dose estimates. One drawback of the PSF method is that it is very difficult to incorporate tissue inhomogeneity into consideration since the point source function can only be generated assuming an infinite, homogenous medium. The point source function is isotropic and solely a function of r, the distance from the point source to target point. Another difficulty is the choice of grid size. In PSF, the doses from source point to target are given as a function of the distance between them. The assumption is that the size of source and target cells are much smaller than the distance. For a radionuclide of interest in RIT, the beta range is usually fairly small (e.g., in the order of mm). Therefore, in order to use the PSF, a patient's CT slice must be subdivided into a small grid matrix so that each pixel is small enough and can be approximated as a point. There are several difficulties in doing this. First, the computation time increases geometrically with the number of points. Sgouros et al. reported an overnight calculation using PSF to estimate dose contours for only one organ (liver) (G. Sgouros, S. Chiu, K. S. Pentlow, et al. "Three-Dimensional Dosimetry for Radioimmunotherapy Treatment Planning," J. Nucl. Med. 34, 1595-1601 (1993)). Second, the nuclear image data acquired usually is much coarser than the calculation grid size. Thus it does not make very much sense to carry the dose calculation to such a fine scale. On the other hand, if the grid size is made too large, e.g., comparable to the size of nuclear medicine image pixels, the error caused by approximating each pixel as a point might be fairly large. This is particularly true when the range of beta emitted is smaller than or comparable to the size of nuclear image pixels. In such a case, the dose to each point is mainly contributed by the cell containing it and close neighboring cells.

As explained above, for radioimmunotherapy, where a beta emitter is used to deliver therapeutic doses to tumor sites, the accurate estimation of the radiation dose to therapy-limiting tissues is much more demanding in order to maximize the therapeutic benefit while minimizing unnecessary radiation exposure.

It is therefore an object of the present invention to provide a system and method that determines internal radioactivity and absorbed dose estimates of a patient more accurately and with greater ease.

It is a further object of the present invention to provide a system and method for determining internal radioactivity of a patient even with high background radioactivity and/or overlapping organs by using uptake information from the gamma camera images and anatomic data from CT scans.

It is another object of the present invention to provide a system and method for determining internal absorbed dose estimates of a patient by using the actual patient's anatomy in terms of organ shape and size as well as inter-organ position, without the use of phantom shapes, as in conventional absorbed dose estimation methods.

It is a further object of the present invention to provide a system and method for determining internal absorbed dose estimates of a patient without assuming a uniform distribution of activity in the source organs, so as to be able to incorporate non-uniform distribution of the source activity if the information is available.

It is a further object of the present invention to provide a system and method for determining internal absorbed dose estimates of a patient which provides not only an average absorbed dose in the target organs but also provides spatially varying radiation dose or statistical information of the dose distribution within an organ (e.g., a volume-dose histogram).

It is a further object of the present invention to provide a system and method for determining internal absorbed dose estimates of a patient in which the dose estimates for beta particles are calculated without assuming local absorption and in which consideration is given to the range of the beta particles; i.e., cross-fire between adjacent organs.

The foregoing and other objects, features and advantages of the invention will be apparent from the following more particular description of a preferred embodiment of the invention, as illustrated in the accompanying drawings wherein like reference numbers represent like parts of the invention.

FIG. 1 is a block diagram of a preferred embodiment of the internal radiation and absorbed dose estimation system of the present invention.

FIG. 2 is a flow chart depicting the steps employed in the present invention to determine the internal radioactivity of a patient.

FIG. 3 shows a flow chart depicting the steps employed in estimating the average dose to an organ and in obtaining statistical information.

FIG. 4 shows a MCNP simulated scale point source dose kernel, F(r/X.sub.90), for .sup.90 Y beta decay, with results of other references shown for comparison.

FIG. 5 shows a volume-dose histogram for liver in a patient.

FIG. 6 shows a volume-dose histogram for kidneys in a patient.

FIG. 7 shows a volume-dose histogram for a right lung in a patient.

FIG. 8 shows whole body nuclear medicine images with regions of interest from the CT data set, from left to right, anterior image, posterior image, and anterior image with regions of interest from CT data.

The present invention may be implemented in a computer system 1 with a CT/MRI system 2 and a gamma camera 3, such as that shown in FIG. 1.

The computer system is preferably a SUN-Sparc 10 workstation, but any workstation, personal computer, minicomputer or mainframe computer may be employed. The computer system preferably has a processing unit, random access memory and permanent storage. Also, an optical communications adaptor may enable communication with other computer systems.

FIG. 1 provides an exemplary arrangement of components used in the preferred embodiment. Therefore, other components with similar function could replace those shown in FIG. 1.

The preferred embodiment is a computer implemented process operating on computer system 1. Images from CT/MRI 2 and Gamma Camera 3 are transmitted to computer 1 to be processed and viewed on monitor 5. A user may operate keyboard 4 as desired. Computer 1 then carries out the operations described below using a computer program product 6 having a computer readable code 7.

The present invention operates according to the steps shown in FIGS. 2 and 3. First, a patient is CT or MRI scanned to obtain CT slices over the relevant regions of the patient's body at 201. For example, a GE Hi-speed Advantage CT system (Milwaukee, Wis.) may be used to acquire CT images. For a patient's chest and abdominal series, a slice interval of 1.0 cm may be used.

Regions of interest are then drawn on the CT slices to obtain a CT volume at 202.

Next, antibodies are labelled with .sup.111 In and .sup.90 Y at 203, and the antibodies are injected into a patient at 204.

Then, an area of the patient is imaged with a gamma camera to obtain anterior and/or posterior gamma images at 205. For example, a Toshiba dual-head scintillation gamma camera (model GCA-7200A) fitted with medium energy parallel-hole collimator may be used to acquire nuclear medicine images. The gamma camera may be peaked for .sup.111 In (photon energies 173 Kev and 245 Kev) with a 20% symmetric window.

Next, the anatomical information from the CT volume and the radioactivity information from the gamma camera image are fused at 206. To do this, the CT volume is transformed into a two dimensional coronal projection and added to the two dimensional gamma camera image. That is, the CT volume is adjusted in size to have the same scale as the nuclear medicine planar images.

More specifically, the CT image is visually moved on a monitor until the organ boundaries on the CT image coincides with those on the gamma camera images. When the organ boundaries are not visible on nuclear medicine images, body contours from both the CT projection and nuclear medicine images could be used as a marker to help registration.

A user controls the fusion operation registering CT and gamma camera images. This is used instead of an autoregistration technique to avoid problems caused by large anatomical changes between CT and gamma camera scans.

That is, during the CT procedure, patients are asked to hold their breath, while patients are allowed to breath normally during the gamma camera scan. This causes the position of the patient's diaphragm to differ between the two images. Therefore, the registration cannot be a simple global translation and rotation operation since the relative inter-organ positions can change between the two imaging modalities.

A second reason for a user controlled fusion approach is that contrast on the gamma camera images is relatively poor so that auto-registration would not work very well, but on the other hand, the visual fusion technique is appropriate for this situation.

Next, to obtain a linear radioactivity density, a set of equations are developed by uniformly selecting several hundred points (N) from the gamma camera image at 207.

More specifically, for each image point (i) selected: ##EQU3## where C.sub.i =cpm from gamma camera image at the calculation point;

a.sub.i (x)=linear radioactivity density at depth of x along the point (.mu.Ci/cm);

B.sub.i (x),.mu. and CF have been defined in equation (1b) above.

Notice that this result is a generalization of equation (1a).

Under the assumption that the radioactivity density, i.e., a.sub.i (x), in any organ is uniform, equation (4) can be simplified by taking the linear radioactivity density out of the integral. It may also be assumed that B.sub.i (x) does not vary appreciably over the range of integration so that it also may be replaced by an average buildup factor. Thus: ##EQU4## where n=the number of organs that can be seen under the calculation point, or in the general case, the number of organs for which we want to calculate the activity. n is on the order of 10 in clinical situations.

a.sub.j =the average linear radioactivity density in the organ j(.mu.Ci/cm). This is the set of unknowns, 1<j<n.

B.sub.j =the average buildup factor for the organ j taken as the buildup factor at the center position of the organ;

l.sub.1j and l.sub.2j =the distances from the upper and lower boundary of the organ j to the body surface respectively (cm).

In equation (5), l.sub.1j and l.sub.2j can be obtained from the CT volume image. So a.sub.j are the only unknown variables (n unknowns). With N points (N>>n) selected, equations (5) form an overdetermined linear equations set.

where ›H! is defined as the product of the known quantities in equation (5). Bold letters are used to represent matrices.

Equation (6) can be solved by a numerical matrix inversion which is equivalent to least square fit of the equation set:

This numerical matrix inversion gives the average activity density for the organs as well as the associated statistical standard deviations for the variables estimated.

Then, to calculate total radioactivity in the organs of interest, the linear activity concentration is multiplied by the CT volume of each organ at 209.

In order to get the effective attenuation coefficient and buildup factor, a thin uniform circular source of diameter 4.89 cm containing 880 .mu.Ci (32.56 MBq) of .sup.111 In may be used. The source is imaged in air and at multiple depths created by tissue-equivalent polystyrene slices, ranging from 1.0 cm to 28.0 cm. The source is made as uniform as possible so that there is a constant amount of radioactivity per unit area.

After the images are taken, the curve representing the buildup factor and calibration factor (CF) of the camera is calculated on a per-pixel basis, i.e., count rate (cpm) per pixel of the image versus radioactivity in each pixel is acquired to be used in the calculation procedures. Using program SAAM II, a buildup factor of the form:

is fit to the data.

The cumulated activity within each organ can be integrated from time-activity data by simple trapezoidal interpolation. A 3-D volume of the patient is built from patient's CT slices. Each voxel in the CT volume may have dimensions of 0.4 cm.times.0.4 cm.times.1.0 cm.

The present invention has been applied to a clinical study using the chimeric anti-CEA monoclonal antibody (cT84.66) with results given in Table 1. The patient (JC) was injected with 5 mCi .sup.111 In labeled Mab. Images shown in FIG. 8 were taken 2.38 hours after injection. For large organs with high activity contrast to background, e.g., liver, the difference between the Thomas' geometric mean method and the CAMI method was generally small. But for other organs with lower activity like kidneys and lungs, the difference could be as large as 50%. For those organs whose boundaries were not visible in the gamma camera images, conventional methods cannot be used to estimate activity. Organs with higher counts, e.g., liver, spleen and heart, usually have smaller standard deviation of the estimated activity due to better statistics. For later time points, the radioactivity left in the organs became less, so the standard deviation became larger.

It can be seen from Table 1 that the CAMI gave internally consistent results. Two parts of the same system, left kidney versus right kidney and left lung versus right lung, have very similar linear density. For example, the left kidney uptake (in FIA) at 2.38 hours was estimated as 0.0123.+-.21.2% while the right kidney value was 0.0102.+-.19.4%. In the corresponding lung comparisons, left lung estimated uptake was 0.0465.+-.8.8% whereas right lung estimated uptake was 0.0530.+-.8.5%. This indicated that the uptake of activity was tissue dependent.

The initial results were obtained using a grid of N points placed uniformly over the nuclear planar image to establish linear densities. No dependence upon N was observed if the number of points was between 50 and 840. Moreover, no improvement in accuracy occurred with the use of a bilateral instead of unilateral nuclear image data. However, background uptake was consistently overestimated via the CAMI methodology in this context of a uniform grid.

A separate experiment was carried out to determine if Compton scatter was enhancing the apparent uptake in the background volume. By positioning a linear (one-dimensional) grid of matrix inversion points at different distances from the edge of a phantom organ with relatively high uptake (spleen), the scatter effect could be reduced by over a factor of two by moving the computational volume some 5 cm or further from the organ edge.

Such edge regions in background activity density computations may be avoided by using a simple algorithm which places an exclusion zone around high activity sources prior to selection of the N point set.

Errors in the CAMI methodology arise from a variety of causes. From the statistical variances seen in the inversion process, it was demonstrated that such uncertainties were generally small compared to the absolute errors seen in the phantom estimates. Because of the scatter effect seen above, the primary experimental problem is probably elimination of scatter radiation from adjacent tissues and background. Since the outline of the organ(s) is simply acquired in the CT scan, exclusion of a scatter zone around a relatively concentrated source of activity should not be difficult. This problem was compounded here since Compton scatter of the higher energy peak of .sup.111 In caused scattered photons to appear within the lower-energy window of the gamma camera spectrometer. Thus, the situation may be simplified by use of a single photon radionuclide and/or an asymmetric energy window on the camera. Some reduction in count rate would ensue in such conditions.

In addition to scatter, there are approximation errors which occur, of necessity in the matrix inversion process. While it is difficult to be precise about their relative magnitude, the two primary contributions to this uncertainty are due to assumptions made in going from equation (4) to equation (5). There, it was assumed that the density a(x) and the buildup factor B(x) could both be taken out of the integral and replaced by mean values. The first of these assumptions follows directly from the lowest-order approximation to a(x); i.e., that a constant activity density may represent a given organ. In the second case, extracting B(x) from the integral is also a lowest-order Taylor's Series approximation to the actual buildup function. In this case, the B(x) function was expanded about its value at the midpoint of the organ. Future improvements to these simplifications to the analytic process could, therefore, involve better representations of these functions within the integral of equation (4).

Geometric data from CT or MRI scans can be used to determine radionuclide uptake in a phantom or patient. Such anatomic information allows planar nuclear images to be interpreted in terms of organ uptakes at depth. These uptakes, determined via the CAMI algorithm, may be in organs which overlap or are not visible due to low contrast in the nuclear projection data. Conjugate view images and auxiliary nuclear imaging are not required. Accuracy of the CAMI method appeared to be within 15% for the phantom organ uptakes independent of projection or number of points (N>50) selected in the matrix inversion process. An additional feature of the method is that the geometric data intrinsic to the patient are immediately available; thus, patient-specific absorbed radiation doses may be estimated without recourse to MIRD or other generic phantoms.

The present invention does not require two conjugate images to be acquired; a single projection is sufficient to determine the unknown linear activity densities. Additionally, the present invention does not require organs to be visible on the nuclear medicine images and, correspondingly, there is no need to draw ROI's on the nuclear images. The CAMI technique is different from the one proposed by Goris et al. where ROIs are still needed to be drawn on gamma camera images and two conjugate images are still required for geometric mean calculation (M. L. Goris, S. A. Knox, K. R. Nielsen, O. Bouillant, "Organ Modeling in the Quantitation of Planar Images for Distribution Studies," Cancer 73(3), 919-922 (1994)).

The internal radioactivity determined by the CAMI technique may then be used to determine absorbed dose estimates as shown in FIG. 3.

In order to get the voxel source kernel (VSK) for .sup.90 Y, an infinite, homogenous soft tissue medium was set up in MCNP at 301. A general-purpose Monte Carlo code, originally developed at Los Alamos National Lab, may be used. Also, the same chemical composition of the soft tissue used in original MIRD phantom for the S value calculation may be used.

Next, the universe is subdivided into a 3-dimensional voxel matrix whose dimensions are the same as the one in the patient's 3-D CT volume (e.g., 0.4 cm.times.0.4 cm.times.1.0 cm) at 302.

The source is uniformly distributed within the central voxel and beta particles are isotropically emitted. The energy of the beta particles is sampled from a .sup.90 Y energy spectrum. Table 2 shows the MAVSK for .sup.90 Y in soft tissue and Table 3 gives the analogous result in lung tissue. Only the data in the first octant is shown because of spatial symmetry. Notice that since the range of betas from .sup.90 Y is only 1.1 cm, dose contribution from the source voxel to distant voxels is approximately zero assuming Bremsstrahlung radiation can be neglected.

Next, assuming unit activity uniformly distributed in the center voxel (position 0,0,0), the total energy deposition to the neighboring voxels as well as the source voxel are simulated in MCNP at 303.

Then, average doses in each voxel are calculated by dividing the energy deposited by the mass of each voxel at 304 and a look-up table is generated (Table 3).

The voxel source kernel for .sup.90 Y in lung tissue (.rho.=0.3 gm/cm.sup.3) is also calculated (Table 3) in a similar manner. In the simulation for the lung voxel source kernel, the chemical composition of lung tissue obtained from any technical reference or source may be used.

Next, a uniform activity concentration in each organ is obtained by dividing cumulated activity by the total number of voxels in that organ at 305. Then, the average absorbed dose to each voxel in the 3-D patient's CT volume is estimated by convoluting the

MAVSK (kernel) with the 3-D activity distribution obtained from CAMI at 306.

Then the average absorbed dose in each organ is calculated by averaging the doses in all the voxels within the organ at 307. Statistical information such as standard deviation of the average dose in each organ as well as the dose-volume histogram is also obtained at 308.

Although uniform activity distribution within each organ is assumed in the above example, a non-uniform distribution of radioactivity could be easily incorporated into the present invention by using a radioactivity concentration look-up table at step 306.

Furthermore, if the tumors are visible in a patient's CT scan, tumor volumes can be obtained by drawing a boundary around the lesion. Thus, absorbed dose to tumor can also be estimated by the present MAVSK system and method.

The present invention was tested by simulating the point source function for .sup.90 Y and comparing the result with several published PSFs. In this example, MCNP was set up in such a way that the electron spectrum of .sup.90 Y is emitted isotopically from a point source in an infinitely homogeneous H.sub.2 O medium. Electrons and secondary photons are followed.

Energy deposited along each electron step is scored in spherical shells of thickness 25 .mu.m. Electrons and photons are transported down to 10 Kev, below which their energy is deposited on-site and the history terminated. This energy cutoff is chosen since the CSDA (continuous slowing down approximation) range of such electrons in H.sub.2 O is exceedingly short (2.5 .mu.m). Two million electron histories are followed, from which the mean energy imparted to each shell and the associated standard deviation are determined. Results are shown in FIG. 4.

The results are shown in the form of scaled absorbed dose versus percentile distance which are defined by Berger. Standard deviations for the MCNP results are all smaller than 1.0% and are not displayed. Data included are from Berger et al., Prestwich et al. and Simpkin et al. (M. J. Berger, "Distribution of Absorbed Dose Around Point Sources of Electrons and Beta Particles in Water and Other Media," MIRD Pamphlet No. 7, J. Nucl. Med., Volume 12: Suppl. No. 5, pp. 5-24 (1971); W. V. Prestwich, J. Nunes, and C. S. Kwok, "Beta Dose Point Kernels for Radionuclides of Potential Use in Radioimmunotherapy," J. Nucl. Med., Volume 30, pp. 1036-1046 (1989); D. J. Simpkin, T. R. Mackie, "EGS4 Monte Carlo Determination of the Beta Dose Kernel in Water," Med. Phys., Volume 17(2), pp. 179-186 (1990)). The simulated results (solid line) from MCNP are in very good agreement with experimental data as well as with other Monte Carlo simulations (e.g., EGS4) with differences at all distances generally within 2%.

A patient study was performed using patient activity data calculated by the CAMI method. The area under the curve (AUC) yielded residence times (in MBqHour) for each organ which could be inserted into the MIRDOSE3 program to obtain dose estimates. All AUCs used in this embodiment are calculated by simple trapezoidal integration. Integration after the last time point is carried out by assuming purely physical decay of the isotope, which is probably an overestimate of the residence time.

Absorbed dose estimates based on AUC data determined above are shown in Table 4. For comparison,

results using the MIRDOSE3 program are also shown. Patient #1 was a normal size female (weight of 114 lbs, height of 5'1), similar to the female phantom described in MIRDOSE3 (nominal weight=125 lbs). Difference between MAVSK and MIRD results were not very marked. Patient #2 was a large size female (weight of 186 lbs, height of 5'6") with some large size organs (e.g., liver and spleen). As can be seen, the differences between MAVSK and MIRDOSE3 estimates could be as large as a factor of three in this case.

A simple modification to MIRDOSE3, suggested by Mike Stabin at Oak Ridge National Lab, is to adjust the MIRDOSE3 results by the patient's organ masses when those organs are very different from standard MIRD sizes. Adjusted MIRDOSE results can be calculated by the following formula; ##EQU5## The results are shown in Table 5. For liver, spleen, kidneys and heart, the radioactivity concentration are higher than neighboring organs and tissue background. The dose contribution are mainly due to the radioactivity inside themselves. Adjusted MIRDOSE results become fairly close to the ones estimated by MAVSK with difference smaller than 15%. However, for lungs, a relatively large portion of dose comes from neighboring organs. The difference is still significant between MIRDOSE3 and MAVSK with errors of about 40%. This occurred because MIRDOSE3 neglects dose contributions from other organs for beta particle irradiations.

An advantage of the MAVSK method over the current MIRDOSE3 program is that it gives statistical information of the absorbed dose distribution. FIGS.

5-7 show the volume-dose histogram (VDH) for patient #2. The present MIRDOSE3 algorithm only gives such information for red marrow absorbed dose. For all the other organs, MIRDOSE3 assumes a Dirac delta function for the VDH; i.e., every voxel in the organ receives a single dose. FIG. 5 shows VDH for the liver. Since the dose to liver voxels mainly came from the activity in liver itself, it was very similar to a delta function distribution. FIG. 6 shows the VDH for both kidneys. As can be seen, the cross-organ dose to kidneys from liver and spleen clearly made VDH farther away from a delta function. Since the right kidney was closer to the liver, some portion of right kidney received higher dose than average. This effect became more marked in the lung case (FIG. 7). Some portions of the right lung (voxels near the liver) had a radiation dose 10 times higher than average absorbed dose. Since the lungs is one of the major dose-limiting organs in RIT, this effect becomes essential during RIT treatment planning. Even if average lung dose is much below the lung dose threshold, some portions of the lungs still get much higher dose which later can become a carcinogenic factor.

A second advantage of the MAVSK method over Monte Carlo Calculation and other point source function methods is its real-time dose estimation capability. The Monte Carlo simulation of the VSK is a separate general calculation which can be run for relatively long periods to get good statistics on the voxel kernels. With a 4 day MCNP run on a SUN-Sparc 10 (Model 51) workstation, standard deviations obtained were less than 2% for all the voxels of interest. After the VSKs were acquired, the patient dose calculation was just the convolution of the VSKs over patient CT volume. This process took less than 1 minute per patient. Thus, patient-specific real-time dose estimates can be accomplished.

For beta radiation, an additional simplifying assumption made in the MIRDOSE3 program is that the beta radiation originating in the source organ is wholly absorbed in that organ without cross irradiation of adjacent tissue. For small or close-neighboring organs and high energy beta particle emitters, this assumption also needs to be improved. This limitation also occurs in its tumor dose estimation method. In the latest release of MIRDOSE3, the algorithm only considers absorbed dose contributed by activity located in the tumor itself. Dose contributions from surrounding organs or even the tumor-bearing organ are not considered. For example, absorbed dose to tumor from liver in the case of hepatic metastases or hepatoma. Similar improvements are also needed in the red marrow absorbed dose estimation to account for radioactivity in adjacent source regions, skeletal structures, and circulating blood. Thus it is clear that the current MIRDOSE system, while not necessarily wrong, needs to be modified so that it can be more flexible and can do better absorbed dose estimation for RIT.

In sum, the above described system and method determines internal radioactivity and absorbed dose estimates with more ease and accuracy than conventional systems.

The CAMI method permits uptake estimation in the common case where two or more organs overlap in the nuclear planar image; e.g., liver and right kidney in a planar projection. It is also applicable in cases where a planar image is simply not obtainable due to a lack of image contrast.

By alignment of nuclear planar and CT images, a simple matrix inversion process determines, via least squares, most likely values for the linear density of activity in the organ(s) of interest. These activity densities were then integrated over the organ volume to determine the total uptake in each tissue.

The Monte Carlo assisted voxel source kernel (MAVSK) system and method provides a novel way to perform real-time patient-specific absorbed doses for beta radionuclides. MAVSK does not require any extra nuclear scans for the patient and is easy to implement. Patient-specific organ doses as well as the dose volume histogram, and even 3-D isodose contours can all be calculated within 1 minute on a SUN-Sparc 10 workstation. Furthermore, a non-uniform distribution of radioactivity can be easily incorporated into MAVSK. Finally, tumor doses can also be estimated if desired; these doses could include contributions from adjacent organs.

The invention having been thus described, it will be apparent to those skilled in the art that the same may be varied in many ways without departing from the spirit and scope of the invention. Any and all such modifications are intended to be included within the scope of the following claims.