Image acquisition and reconstruction
The time course of a radioactive concentration in a volume was sampled at multiple time points, using a Siemens ECAT HR+ (PETonly) scanner [17], a General Electric (GE) Discovery ST (PET/CT) scanner [18], or a GE Discovery 690 (PET/CT) scanner [19]. The described methods have been employed in a total of about 100 scans, with a large number of different tracers and diseases. A subset of these has been analyzed for this paper. The acquisition data for the cases reported in this paper will be described below.
Acquisition 1: Image data from one patient with a known cancer in the headneck area was used. This study was part of a larger project in which we sought to evaluate the effect of treatment on cancer metabolism using serial [^{11}C]acetate PET scans. As part of this project each PET study consisted of a short dynamic scan over the heart region, immediately followed by a longer dynamic scan over the area of cancerous growth. The heart scan was used to establish an imagederived input function for subsequent kinetic modeling of acetate uptake in the tumors. Cardiac images were acquired on the GE Discovery ST [18], using 50 MBq [^{11}C]Acetate, which is only 5% of the standard dose for [^{11}C]acetate cardiac scanning in our institution. Acquisition comprised a 23 frame sequence (12×5s, 6×10s, 4×30s, 1×60s) started at injection. Images were reconstructed to a 30 cm display field of view employing standard Filtered Back Projection (FBP), with a 6.3 mm Hanning post filter.
Acquisition 2: Following the heart scan, headandneck images were acquired on the GE Discovery ST [18], using 1000 MBq [^{11}C]Acetate, with a 32 minute frame sequence (12×5s, 6×10s, 4×30s, 4×60, 2×120, 4×300) starting at injection. Images were reconstructed to 30 cm display field of view, using the above FBP algorithm. The same data was also reconstructed to 30 cm display field of view, using an ordered subsetsexpectation maximization (OSEM) algorithm (21 subsets, 2 iterations, 3.91 mm FWHM loop filter, and 5.45 mm postfilter).
Acquisition 3: The brain images (used as one of the tracers comparing 2D and 3D Hotelling filter) were acquired on healthy volunteers, on the GE Discovery 690. A total of 250 MBq [^{11}C]Raclopride [20] was injected as bolus+infusion [21]. The 21 frame scan (9 frames × 2 min , 3 frames × 3 min, 3 frames × 4.20 min, 3 frames × 5 min) was acquired from start of injection, employing the GE Discovery 690. Images were reconstructed to a 30 cm display field of view employing the iterative VuePoint HD algorithm, with 2 iterations 24 subsets and 6.4 mm postfilter.
Acquisition 4: The brain images (used for Patlak parametric images) were acquired on the Siemens ECAT HR+ scanner [
17], using 200 MBq FDG. A 90frame scan (90 × 30s) was acquired starting at injection. Images were reconstructed to a 30 cm display field of view employing FBP, with a 4 mm Hanning post filter. The 90 images were grouped sequentially in 9 groups each containing 10 images. New 9frame dynamic sequences were created from these groups by:

Selecting one image from each group, creating a 9frame dynamic scan with statistics mimicking 10% (20 MBq) of the injected activity.

Summing three images from each group, creating a 9frame dynamic scan with statistics mimicking 30% (60 MBq) of the normally injected activity

Summing all images in each group, creating a 9frame dynamic scan from the 90frame scan with statistics equal to 100% (200 MBq), that is, the complete injected activity.

The plasma input function was measured with arterial blood sampling.
Acquisition 5: The liver images (used for Patlak parametric images) were acquired on the GE Discovery ST using 370 MBq [^{18}F]FDG in a patient with cholangiocarcinoma and liver metastases. The scan included 16 frames (5× 60s, 5×180s, 6×300s ), acquired during 50 minutes employing 2D mode on a General Electric Discovery ST (PET/CT) scanner [18]. Images were reconstructed to a 50 cm display field of view employing FBP, with a 4 mm Hanning post filter. The plasma input function was measured with arterial blood sampling. Also an alternative image derived input function was calculated using a manually delineated regionofinterest in multiple slices of the ascending aorta, which was visible in the first frames, and clearly visible when applying the Hotelling filter.
Acquisition 6: Dynamic imaging of about 200 MBq [^{18}F]fluorothymidine (FLT) was used for determination of typical image noise levels, employing a GE Discovery 690 [19] in timeofflight mode.
The image material was obtained from patient studies in research projects that were conducted in accordance with the declaration of Helsinki and approved by the local human ethics committee and the hospital radiation ethics committe (UppsalaÖrebro EPN reference 2006/130932, Uppsala University Hospital reference SEK 2005:03, and Umeå University Hospital reference 201126331M). All subjects provided written consent prior to the study.
Image processing
The new filter uses the PCA transform [16], also known as the Hotelling transform. The Hotelling filter can be applied to either the data set containing all time frames, or a subset where early frames are excluded from the principal component transforms, as described below. In the case when a subset is filtered, the final resulting data set is assembled by concatenating the original excluded timeframes with the Hotelling filtered timeframes. Thus, the resulting filtered data set will also in this circumstance contain the same number of frames as the original data set. We will now give a detailed description of the Hotelling filtering algorithm:
1) From image to principal component domain:
We organize the dynamic data for a single tomographic slice (2D Hotelling) or a volumetric matrix (3D Hotelling) into a 2dimensional matrix, with the row vectors being pixel values for each of the
N timeframes. We exclude the pixels that have zeros in all frames (columns with data points being all zero). Inspired by [
16], each row (image) in this matrix is standardized by subtracting its mean and dividing with its sample standard deviation, giving the matrix
X. A new matrix
A is formed by sorting the eigenvectors of the covariance matrix of
X, in descending order according to the corresponding normalized eigenvalues
λ
_{
k
} (one per eigenvector). The PCA transform [
16] gives weighted sums
Y of the original dynamic data, called “principal components”:
As an option, at this point “principal component” images can be created. This is done by putting back the excluded columns (pixels that had zeros in all frames) in their correct position in the Y matrix, and reorganizing the rows to twodimensional images.
2) Remove noise:
The eigenvalues describe the variance of the data which is explained by the principalcomponent images (stored as rows in) Y. The eigenvectors with low eigenvalues are removed from A, setting the corresponding vector values to zero, giving the matrix B.
Thus, the first principal components 1 to n are used, and the components n+1 to N are set to zero. The Hotelling filter will be described as PC1n (for instance PC14 if the first n=4 components are used).
3) From principal component domain to image domain:
The PCA transform is reversed to give an approximation
X’ of the original images
using that B
^{1}=B
^{T} because B is a matrix of orthonormal eigenvectors [15].
The standardization of the pixel values are inversed, by multiplying the rows in X’ by the previously calculated values of the sample standard deviation and adding the sample mean. The pixels excluded from the data vector in step 1 above (pixels that were zero in all frames) are put back in their original position. Finally, the pixels of each frame are reorganized to form a set of filtered twodimensional images (2D Hotelling) or threedimensional (3D Hotelling) images, where the number of filtered images is the same as the number of original timeframes.
The eigenvectors and eigenvalues were calculated using the “eig” function in Matlab, and the vectors were sorted according to their corresponding eigenvalues.
For a 27frame, 128*128 matrix dynamic scan, producing a Hotelling filtered image takes 0.08 seconds, on a 2.6 GHz Pentium Dualcore E5300 processor.
An explanation factor
e
_{
k
} that describes the percent of the original variance, accounted for in the
k:th principal component, is calculated
The “explained variance” is defined as the sum of the explanation factors over the n used principal components.
explained variance=
The explained variance is thus a measure of how much of the original variance in the dynamic image sequence remains after filtration. The suggested interpretation is that filtering two image sequences using the same number of principal components, the explained variance will be higher for the image sequence containing less noise. For a dynamic imagesequence, the explained variance increases when increasing the number of principal components employed in the Hotelling filter. The explained variance for the original image (employing all components) is 100%.
The explanation factor and inspection of the principal component images was used to select the correct number of principal components. The number of principal components to include in the Hotelling filter was determined using a Scree plot [22], that is plotting the explanation factor (or normalized eigenvalue) as a function of highest principal component (PC) number employed in the Hotelling filter. The explanation factor decreases with increasing PC number, and abruptly converges to a level close to zero. We have found that taking the PC number where the explanation factor converges to zero, and adding one gives good values for the Hotelling filter.
We have also created a second type of Scree plots, by plotting a “ROI explanation factor”, using the pixel values in the principal component images. Thus for each principal component image and ROI, a ROI explanation factor was formed as the average of the pixel values within a ROI. A Scree plot was formed by plotting these ROI explanation factors as a function of component number.
As a quality control, the residual image, defined as the difference between original image and filtered image, was inspected.
Regionof interests (ROI) were analyzed in the residual images creating plots that describe the removed uptake in a ROI as a function of time. The purpose of that kind of analysis is to find bias as a function of time. The number of principal components employed in the comparison between the 2D and 3D Hotelling filters were PC14.
The filtering software and Patlak slopes software was implemented as part of an inhouse developed pixelbased 4dimensional VOI tool, called imlook4d 2.00 [23], written in MATLAB 7.0 (The Mathworks, Inc., Natick, MA). The Hotelling filter is implemented as an integral part of imlook4d, giving the operator interactive feedback to any change in filtering parameters. Similarly, the principal component images and the residual (that is, the removed noise) can be viewed interactively. The parametric images were created with imlook4d, also here recalculated “onthefly”, to get an interactive feel when varying the 2D Hotelling filter. The 3D Hotelling filter was applied using a script to calculate a new filtered matrix.
Patlakslope images were calculated using data starting at 20 minutes past injection. It was checked that data past 20 minutes gave a linear Patlak plot.
The displayed images and small image inserts are always displayed with equal color scales, if not explicitly stated otherwise.
Simulations
A simulation study was performed to study dynamic data that follows a relatively large number of rate constants, a scenario we could not create with a simple radioactive phantom. We have chosen to simulate the uptake by describing typical blood and tissue uptake by mathematical functions describing the blood peak shape, and an uptake curveshape similar to a capacitor charging curve. The reason that we did not create these by compartment modeling is that we did not have access to a true blood input functions for different tracers with varying uptake rates. Although this method can’t be claimed to be an exact representation of a specific tracer, we believe that the method allows us to freely vary time constants and thus in a simple way present more complex data than is typically present in a PET scan. Furthermore, we chose not to mimic any realistic geometry, since we did not want the reader to believe that the simulation is specific to any particular tracer uptake. Also, there is nothing in the mathematics behind the Hotelling filter that correlates with the pixel location, so we prefer to make square regions. Additional file 1 displays an example of the simulated data, and the generated activity curves.
The procedure to create this data with kinetic and noise properties similar to a dynamic PET scan will now be described. A dynamic image sequence was created consisting of 20 noiseless images, each 90 by 90 pixels. In these images, analogous to timeframes in a real PET scan, nine square areas were created, each 30 times 30 pixels in size. Each area was assigned a kinetic behavior similar to blood or tissue uptake. The tissue uptake was simulated to give different kinetic risetimes and uptake amplitudes, using the formula
with permutations of the values A=1,2,3, and k=1,2,3, where q is a factor that was calculated as described below. The image number i is an integer between 0 and 19, which represents time.
The blood uptake was calculated using the formula
The number of blood areas differed in different simulations, since we wanted to study the impact of the blood signal. We employed k=1 for the first blood area, k=2 for a second blood area, and k=3 for a third blood area, when applicable. The factor B was calculated as described below. The number of blood areas was varied between zero and three, replacing tissue areas in the image as necessary to get nine areas in the image.
Gaussiandistributed noise with a standard deviation
was added to the pixel values, with X
_{
i,j
} being the noiseless pixel values in image i and square area j. The noise was created using the Matlab “randn” function.
The factors q and B that gave the most PETlike images were determined as follows. Regionofinterests (ROIs) were drawn over tissue, lung, and blood in two FLT and two FDG dynamic scans. The mean pixel value
and the standard deviation s
_{
i,j
} was calculated for each timeframe image i and ROI j. The relative standard deviation
was found to consistently be in the range between (0.18±0.05) to (0.58±0.10) for all four scans, calculated for each individual timeframe and ROI. The factors q and B were optimized in integer steps so that the simulated dynamic sequence matched this range of relative standard deviations. The values q=10 and B=5 gave relative standard deviations between 0.17 and 0.64 calculated for each individual test image area, and image.
Analyzing the consequences of Hotelling filtering on simulated data was done by applying a circular ROI consisting of 305 pixels within each simulated square area. The reason for using circular ROIs was only practical, since that is the default ROIshape in imlook4d. It is plausible that filtering introduces a bias varying between early and late images, and between ROIs representing different kinetic behavior. Therefore a measure was introduced to describe the timevarying difference between noisy data (filtered or not filtered) with that of the noiseless data. This measure should have the properties that a negative and a positive difference can not cancel, so we calculated the absolute value of the noisy ROI value
x
_{
i,j
} and noiseless
X
_{
i,j
} ROI value for each ROI
j and frame
i. In order to get one value that summarized the bias for all different kinetics and summed over all frames we introduced
with N=180 (20 images times 9 ROIs).
The amount of noise in Hotelling filtered data was quantified as the standard deviation
s
_{
i,j
} of pixel values within ROI
j from image
i. The average standard deviation for all images and ROIs was calculated as
The amount of noise in the orginal noisy data was calculated using the same formula on the original data, but is referred to as
S. The ratios of the average standard deviations of filtered and original noisy data describe the amount of noise that remains following the filtering procedure:
The spatial resolution was quantified by measuring the fullwidth half maximum (FWHM) under different conditions. The simulated images (which have step function interfaces between the square areas) were given a finite resolution by convoluting the image with a Gaussian shaped 2dimensional kernel with FWHM of 4 pixels. After applying the Gaussian filter, the steplike intensity profile along a line crossing an interface thus becomes a smoothed transition between the intensity levels of the squares. Differentiating this profile gives a Gaussianshaped profile. The FWHM of the differentiated profile was measured manually in a graph, and compared for different combinations of Hotellingfilter settings. Average values for the resolution were calculated, with error estimates being the standard deviation of these measures.