PET image reconstruction with rotationally symmetric polygonal pixel grid based highly compressible system matrix

PET image reconstruction with rotationally symmetric polygonal pixel grid based highly compressible system matrix 1 first-author 1 2 2 1 1 1 corresp Corresponding author. E-mail address: maty@tsinghua.edu.cn Corresponding author. E-mail address: maty@tsinghua.edu.cn maty@tsinghua.edu.cn Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry of Education Department of Engineering Physics, Tsinghua University, Beijing 10084, China Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry of Education Department of Engineering Physics, Tsinghua University, Beijing 10084, China National Institute of Metrology, Beijing 100013, China National Institute of Metrology, Beijing 100013, China To achieve a maximum compression of system matrix in positron emission tomography (PET) image reconstruction, we proposed a polygonal image pixel division strategy in accordance with rotationally symmetric PET geometry. Geometrical definition and indexing rule for polygonal pixels were established. Image conversion from polygonal pixel structure to conventional rectangular pixel structure was implemented using a conversion matrix. A set of test images were analytically defined in polygonal pixel structure, converted to conventional rectangular pixel based images, and correctly displayed which verified the correctness of the image definition, conversion description and conversion of polygonal pixel structure. A compressed system matrix for PET image recon was generated by tap model and tested by forward-projecting three different distributions of radioactive sources to the sinogram domain and comparing them with theoretical predictions. On a practical small animal PET scanner, a compress ratio of 12.6:1 of the system matrix size was achieved with the polygonal pixel structure, comparing with the conventional rectangular pixel based tap-mode one. OS-EM iterative image reconstruction algorithms with the polygonal and conventional Cartesian pixel grid were developed. A hot rod phantom was detected and reconstructed based on these two grids with reasonable time cost. Image resolution of reconstructed images was both 1.35 mm. We conclude that it is feasible to reconstruct and display images in a polygonal image pixel structure based on a compressed system matrix in PET image reconstruction. Positron emission tomography (PET) Rotationally symmetric polygonal pixel grid Compressed system matrix Image reconstruction 1 Introduction 1 Nowadays, as an essential branch of molecular imaging technologies, positron emission tomography (PET) has made rapid developments. Based on modern high-sensitivity 3D data acquisition mode and statistical iterative reconstruction algorithms, high-quality molecular and functional images from organs to cells of human bodies can be obtained. PET plays an important role in oncology, cardiology and neurology diagnosis, including early detection, stage of malignant evaluation, and post-treatment monitoring of human diseases. It remains a great challenge though to implement a high-resolution, fully 3D PET reconstruction algorithm in a clinical PET system, since there are numerous number (about 10 7 ) [ 1 1 ] of lines-of-response (LORs) events need to be calculated when a modern PET machine works in 3D mode. Correspondingly, a huge-sized (about 10 12 elements) system matrix is needed in image reconstruction. In order to address this issue, it is critical to compress the system matrix and obtain best trade-off between image quality and reconstruction speed. Considering that most modern clinical PET systems employ a symmetric ring-shape geometry, a concept of taking advantage of this symmetry with irregular image pixel shape instead of conventional rectangular pixels to achieve a maximum compression of system matrix was put forward, which recently becomes a new focus area in fully 3D PET iterative algorithm development. Compared with other methods of system matrix compression [ 2 2 ] , compression with symmetric is direct and has been proven to be efficient to achieve better compression factor. Early approaches designed a polar grid based system matrix structure which has rotational and reflectional symmetries. For example, such principle was proposed to accelerate PET iteration process by polar pixel structure [ 3 3 , 4 4 ] . Recent years, with the popularization of 3D PET technology and explosive growth of acquisition data size, more attention was involved in this field. In addition to the polar grid, various irregular pixel shapes were presented in commercial or laboratory PET systems such as in MicroPET Focus 120 [ 5 5 ] and in LabPET [ 6 6 ] . A software toolkit PRESTO [ 7 7 , 8 8 , 9 9 , 10 10 , 11 11 ] was developed, which integrated the generic calculating, storing and accessing techniques based on the rotational symmetry in a C++ library to improve the efficiency and accuracy of image reconstruction. Comparing to the original polar pixel based geometric structure, an irregular division can overcome the unevenness of pixel area distribution across the image FOV. It has been demonstrated that a polar voxel with square-like shape produces fewer artifacts [ 12 12 , 13 13 ] , so one usually has to design specific division mode for each compression factor as shown in Ref.[ 11 11 ]. The increased complicity of pixel shapes in the irregular division mode brings further challenges in image reconstruction and display. Owing to that most existing PET systems have a polygon ring shape rather than a circle, a polygonal pixel grid may have an edge than traditional polar grid. This work proposes a polygonal image pixel division strategy based on the characteristics of rotationally symmetric PET geometry. It is adaptable to the structure symmetries of various types of PET systems. One can simplify the calculation of the system matrix at a compressed format by the usage of rotational symmetry of pixel structure during the PET image formation process. The proposed method was applied to a practical small animal-PET system. In view of the definition of this pixel grid and calculation of the system matrix, we accomplished the iterative reconstruction and image display of experiment data of a hot rod phantom with 30 mm both in diameter and length based on the practical PET scanner. 3D detected data were organized into sinograms. System matrix was calculated by tap model and compressed by rotational and axial symmetries. A Cartesian grid based image reconstruction algorithms with same pixel size was also generated with tap-model. Storage size of system matrix, time cost and image quality were compared for these two methods. 2 Methods 1 2.1 Rotationally symmetric polygonal pixel based image expression & display 2 The core concept of this paper is based on a novel rotationally symmetric polygonal pixel grid. A representative division pattern of the polygonal pixel grid is shown in Fig.1 Fig.1 (a). The structure contains 8 segments in Fig.1 Fig.1 (a) which are rotationally symmetric with respect to the center of the image. The arrangement of polygonal pixel is same from segment to segment. Compared to the conventional rectangle pixel division pattern ( Fig.1 Fig.1 (b)), the introduction of rotational symmetry helps to compress the storage space of the system matrix. Compared to the existing polar pixel grid pattern ( Fig.1 Fig.1 (c)), the regularity of this division mode ensures a simplified expression and automatic segmentation. In another aspect, there are only three types of basic pixel shapes. One can easily adjust the area of each type of pixel to minimize the difference of pixel area over the image FOV. Therefore, this structure may have an edge over polar grid in terms of reducing pixel area unevenness. Enslaved by the typical image display mode of modern computer, the main challenge of this novel structure is that the image expressed in this rotational symmetry polygon pixel grid has to be translated into the cubic ones for image display. Besides, a new method is needed to calculate geometrical weights of system matrix since the universal algorithms such as Siddons ray-tracing method [ 14 14 ] is not applicable in this particular case. 2.2 Geometrical description of polygonal pixels 2 Based on the concept of building arotationally symmetric geometrical structure, the definition of polygonal pixels consists of three steps: (1) segment definition, (2) layer definition inside segments and (3) pixel definition, as illustrated in Figs.2(a), (b) and (c), respectively. Segment definition: When the compression factor is P , P rays are drawn starting from the center of image field of view (FOV), which divides the FOV into fan-shaped segments. The angle θ between the two adjacent rays, i.e. the vertex angle of each segment is θ =2π/ P . The shape and area of every arrangement is equal to each other. Starting from the positive x -axis which is defined to be the starting boundary of the first segment, an index number p =0, 1… P -1 is anticlockwise assigned to each of the P segments. The ending boundary of p -th segment is the same as the starting boundary of ( p +1)-th segment, where p =0, 1… P -2. The ending boundary of segment P -1 coincides with the positive x -axis. As all the segments are identical in shape, the arrangement of pixels can be defined only once for one segment, which can be rotationally replicated to other segments. Therefore, the definition of layers and pixels are discussed only in segment 0. Layer definition: The area inside one segment is divided into several regions by a group of layer-boundary lines which are perpendicular to the vertex angle bisector and intersect the segment-boundary with equal distance ρ . The distance between two adjacent layer-boundary lines is <math id="M4"><mrow><mi>ρ</mi><mi>cos</mi><mo stretchy="false">(</mo><mi>θ</mi><mo>/</mo><mn>2</mn><mo stretchy="false">)</mo></mrow></math> . The area embraced by the two adjacent parallel lines and the two segment-boundary rays is called a layer, except the centermost layer which is formed by the two segment-boundary rays and one layer-boundary line which is the closest to the vertex angle. Each layer is identified by an index number s . The layer containing the vertex is layer 0, whose adjacent layer is layer 1, and so on. The total number of layers is defined as S , which is determined by the size of the image. 1 1 Parameters 1 1 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10356991&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10356997&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10356992&type= 1 1 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10357001&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10357009&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10357005&type= 1 1 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10357013&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10357020&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10357017&type= 1 1 Shape name 1 1 Triangle 1 1 Right angle trapezoidal 1 1 Rectangle 1 1 Number of pixels 1 1 1 1 1 2 S 1 1 ( S -1) S /2 1 1 Area 1 1 ρ 2 sin θ /2 1 1 3 ρ 2 sin θ /4 1 1 ρ 2 sin θ 1 1 Relative area ratio 1 1 2 1 1 3 1 1 4 Pixel definition: As illustrated in Fig.2(c), a set of lines parallel to the vertex angle bisector are drawn which divide each layer into one or multiple polygonal regions. Those lines are denoted as lay-intersecting lines to be distinguished from layer-boundary line in layer definition, and each polygonal region is defined as one pixel. The distances between the neighborhood layer-intersecting lines are all equal to 2 ρ sin θ /2. The s -th layer has s+ 1 pixels that are in a row. The row-index for those pixels labeled as t . We define that in each layer, the 0 -th pixel contains the ending boundary, and the s -th pixel contains the starting boundary. With the above steps described, one can establish a rotationally symmetric polygonal pixel grid and enumerate all the pixels with three parameters, i.e. segment-index p , layer-index s and row-index t . Therefore, a pixel is uniquely specified by ( p,s,t ). In comparison with conventional rectangular-pixel based image division pattern, the proposed pixel division strategy has a P -fold rotational symmetry when there are P segments, leading to a 1/ P storage size of the system matrix. The existing polar-pixel based methods can achieve a similar compression factor. However, considerable difference of both pixel shape and area exists in pixel division pattern. E.g. in Ref.[ 15 15 ], the largest area ratio of pixels reaches almost 5, which may lead to non-uniform image quality distribution and even image artifacts. As listed in Table 1 Table 1 , in the proposed polygon pixel division strategy, there are only 3 basic pixel types: triangle, right angle trapezoidal and rectangle, and the area ratios of these three types of pixel are 2:3:4. The simple pixel shapes and reduced pixel area difference are beneficial in minimizing resolution non-uniformity and reducing image artifacts. 2.3 Pixel indexing and addressing 2 In order to store a digital image into computer memory, a pixel indexing rule is needed to map each pixel to a memory unit. Likewise, a pixel addressing method is needed to calculate the pixel coordinate given its memory address. In an image expressed by rotationally symmetric polygonal pixels, each pixel is uniquely represented by its three indices ( p,s,t ). When the image is stored in computer, it occupies a continuous space of computer memory, for each pixel ( p,s,t ), a unique address i needs to be determined. Figure 4 Figure 4 shows a representative scheme for assigning memory address to image pixels in an 8-segment, 5-layer image. The pixel indices are determined with a segment-by-segment and layer-by-layer sequence. The pixel indices inside p -th segment are from 15 p to 15 p +14. Inside a segment, the pixel indices are increasingly assigned from smaller layers to larger layers and from smaller rows to larger rows inside each layer. According to above rules, one has the mapping relationship between i and ( p,s,t ) is: <math display="block" id="M5"><mrow><mi>i</mi><mo>=</mo><mi>I</mi><mo>×</mo><mi>p</mi><mo>+</mo><mfrac><mrow><mi>s</mi><mo stretchy="false">(</mo><mi>s</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mn>2</mn></mfrac><mo>+</mo><mi>t</mi></mrow></math> Where <math id="M6"><mrow><mi>I</mi><mo>=</mo><mo stretchy="false">(</mo><mi>S</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mi>S</mi><mo>+</mo><mn>2</mn><mo stretchy="false">)</mo><mo>/</mo><mn>2</mn></mrow></math> is the total number of pixels in one segment. Meanwhile, the coordinate ( p,s,t ) of the i -th pixel is determined as: <math display="block" id="M7"><mrow><mi>p</mi><mo>=</mo><mrow><mo>[</mo><mrow><mfrac><mi>i</mi><mi>I</mi></mfrac></mrow><mo>]</mo></mrow></mrow></math> <math display="block" id="M8"><mrow><mi>s</mi><mo>=</mo><msub><mi>s</mi><mi>t</mi></msub><mo>,</mo><mfrac><mrow><mi>s</mi><mo stretchy="false">(</mo><mi>s</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mn>2</mn></mfrac><mo>≤</mo><mi>i</mi><mo>−</mo><mi>I</mi><mo>×</mo><mi>p</mi><mo>≤</mo><mfrac><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mi>s</mi><mo>+</mo><mn>2</mn><mo stretchy="false">)</mo></mrow><mn>2</mn></mfrac></mrow></math> And <math display="block" id="M9"><mrow><mi>t</mi><mo>=</mo><mi>i</mi><mo>−</mo><mi>p</mi><mo>×</mo><mi>I</mi><mo>−</mo><mfrac><mrow><mi>s</mi><mo stretchy="false">(</mo><mi>s</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mn>2</mn></mfrac></mrow></math> 2.4 Image pixel structure conversion 2 To view an image stored in the rotationally symmetric polygonal pixel structure on computer screen, it has to be firstly converted to a conventional rectangular pixel structure. The conversion process can be expressed as: <math display="block" id="M10"><mrow><msub><mrow><mtext>Image</mtext></mrow><mi>M</mi></msub><msub><mrow><mtext>=Q</mtext></mrow><mrow><mi>M</mi><mi>N</mi></mrow></msub><msub><mrow><mtext>Image</mtext></mrow><mi>N</mi></msub></mrow></math> Where Image M is the converted image in rectangular pixel based structure, Image N is the original image in rotationally symmetric polygonal pixel structure, both organized in column vector format, and Q MN is the conversion matrix. As illustrated in Fig.5 Fig.5 , the value of m -th rectangular pixel in Image M can be expressed as a weighted sum of the values of multiple polygonal pixels overlapping it. The overlap area of the n- th pixel in polygonal structure and the m- th pixel in rectangular structure is denoted as S mn . The weight is considered as the ratio of the intersection area and the total area of the cubic pixel. Therefore, one has <math display="block" id="M11"><mrow><msub><mi>b</mi><mi>m</mi></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>N</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>S</mi><mrow><mi>m</mi><mi>n</mi></mrow></msub><mo>×</mo><msub><mi>a</mi><mi>n</mi></msub></mrow></mstyle></mrow><mrow><msub><mi>S</mi><mi>m</mi></msub></mrow></mfrac><mi>，</mi><mtext> </mtext><msub><mi>q</mi><mrow><mi>m</mi><mi>n</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>S</mi><mrow><mi>m</mi><mi>n</mi></mrow></msub></mrow><mrow><msub><mi>S</mi><mi>m</mi></msub></mrow></mfrac></mrow></math> 2.5 Verification of image expression and structure conversion 2 In order to verify that an image is correctly stored in a rotationally symmetric polygonal structure and converted to a rectangular pixel structure, a set of test images are analytically defined: <math display="block" id="M12"><mtable columnalign="left"><mtr><mtd><mi>a</mi><mrow><mo>(</mo><mrow><mi>p</mi><mo>,</mo><mi>s</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>,</mo><mtext> where </mtext><mi>s</mi><mtext> is odd</mtext></mtd></mtr><mtr><mtd><mo>=</mo><mn>1</mn><mo>,</mo><mtext> where </mtext><mi>s</mi><mtext> is even</mtext></mtd></mtr></mtable></math> By selecting different segment number, a set of images A n is numerically defined. Fig.6 Fig.6 illustrates the test images of 8, 9 and 12 segments, respectively. Those images are converted to B m , stored in rectangular pixel structure, and displayed on computer screen for examination. If the structure is defined correctly and the transforming matrix works effectively, one should see the images corresponding to the predefined ones in various rotationally symmetric polygonal grids. 3 Generation and verification of system matrix 1 3.1 Generation of a compressed system matrix 2 In the calculation of the system matrix, each line-of-response (LOR) is considered as a tape with certain length. As shown in Fig.7 Fig.7 , each LOR is uniquely identified with its projection angle φ and the distance to the center of FOV r . Therefore, the j -th LOR is denoted as j ( φ , r ). Each element c i,j in the system matrix records the probability that the j -th LOR detects a photon emitted from the i -th pixel, denoted as i ( p , s , t ). In this study, the impact of physical effects, such as positron range, photon non-collinearity, attenuation and scattering are ignored, and pure geometric detector response is modeled. In this way, c i,j is considered to be proportional to the overlapping area of the i -th pixel and j -th LOR. It is clear the rotational symmetric structure in the image domain can be utilized to compress the system matrix. As shown in Fig.8 Fig.8 , considering two image pixels i 0 and i , when they have same layer index s and row index t , and when the angle from LOR j 0 to j is equal to pθ , i.e. the angle from segments p 0 to p , the overlapped area between pixel i and LOR j is rotationally symmetric to that between pixel i 0 and LOR j 0 , thus the corresponding system matrix element has same values, i.e. <math id="M13"><mrow><msub><mi>c</mi><mrow><msub><mrow/><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msub></mrow></msub><mo>=</mo><msub><mi>c</mi><mrow><msub><mi>i</mi><mn>0</mn></msub><mo>,</mo><msub><mi>j</mi><mn>0</mn></msub></mrow></msub></mrow></math> . The relationship between i and i 0 is shown followed: <math display="block" id="M14"><mrow><msub><mi>i</mi><mn>0</mn></msub><mo>=</mo><mi>i</mi><mo>−</mo><mi>I</mi><mo>×</mo><mi>p</mi></mrow></math> The relationship between j and j 0 can be expressed by relationships of r and r 0 as well as φ and φ 0 : <math display="block" id="M15"><mrow><msub><mi>r</mi><mn>0</mn></msub><mo>=</mo><mi>r</mi></mrow></math> <math display="block" id="M16"><mrow><msub><mi>φ</mi><mn>0</mn></msub><mo>=</mo><mrow><mo>{</mo><mrow><msubsup><mrow/><mrow><mi>φ</mi><mo>−</mo><mi>p</mi><mi>θ</mi><mo>+</mo><mi>π</mi><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>φ</mi><mo>−</mo><mi>p</mi><mi>θ</mi><mtext> </mtext><mo><</mo><mn>0</mn></mrow><mrow><mi>φ</mi><mo>−</mo><mi>p</mi><mi>θ</mi><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>φ</mi><mo>−</mo><mi>p</mi><mi>θ</mi><mtext> </mtext><mo>></mo><mn>0</mn></mrow></msubsup></mrow></mrow></mrow></math> 3.2 Verification of system matrix 2 The following verification studies are based on a simulated typical small animal PET system. As shown in Fig.9 Fig.9 , the detector consists of 6 rings, on each of which 32 detector blocks are evenly distributed. Every block is composed of 8×8 LYSO crystals. The size of crystal is 1.59 mm×1.59 mm×10 mm. The detector ring diameter is 147.2 mm, and the size of FOV is 30 mm in diameter and 76.32 mm in axial length. Radial center of the entire detector field of view is a circular with a radius of about 30 mm. A rotationally symmetric polygonal pixel structure with P =12, ρ =0.5 mm is defined. The layer-height ρ is chosen so that it is much smaller than the desired PET system resolution i.e. 1.6 mm. To cover the 30 mm FOV, the image has 60 layers. The acquired 2D PET data are organized to a sinogram with 155 angular samplings over 180 degree and 60 radial samplings over 30 mm trans-axial FOV, and the width of the LOR tape is the same as the layer-height. By defining different source distributions in the image structure, the corresponding sinogram data were generated by projecting the source image into LOR space using the pre-calculated system matrix. The sinograms are visually examined and compared to theoretical expectation to verify the correctness of the system matrix and projection model. 4 Image reconstruction 1 Data for image reconstruction were acquired with an experimental small animal PET system. The detector consists of 5 rings, on each of which 24 detector blocks are evenly distributed. Every block is composed of 9×9 LYSO crystals. The size of crystal is 2.11 mm×2.11 mm×10 mm. The inner detector ring diameter is 76.33 mm, and the size of FOV is 50 mm in diameter and 95 mm in axial length. A rotationally symmetric polygonal pixel structure with P =12, ρ =0.33 mm was defined. The layer-height ρ was chosen so that it was much smaller than the desired PET system resolution. To cover the 50 mm FOV, the image had 79 layers. The acquired 3D PET data were organized to a sinogram with 108 angular samplings over 180 degree and 70 radial samplings over 46.76 mm trans-axial FOV, and the width of the LOR tape was the same as the layer-height. After the establishment of the system matrix, the projection data acquired from the real PET scanner were reconstructed based on the rotationally symmetrical structure. Using the image transformation matrix defined in Eq.(5), the reconstructed images were transformed and displayed for visual examination. We also conducted the image reconstruction based on a Cartesian pixel grid with 142×142 pixels of 0.33 mm×0.33 mm area. The system matrix was calculated using the same tap model for LORs. The commonly used OS-EM algorithm [ 16 16 , 17 17 ] was chosen for PET image reconstruction in this study. The OS-EM algorithm is based on maximum likelihood estimation of Poisson distributed signal. The objective of the iterative algorithm is to maximize cost function which defines a logarithmic likelihood function between the projection of the image estimate and the actual projection values. An iterative expectation maximization algorithm is used to optimize the cost function. The projection views are grouped in L subsets, the image is updated after each subset is considered. The iterative formula is as followed: <math display="block" id="M20"><mrow><msubsup><mi>f</mi><mrow><mi>s</mi><mi>i</mi></mrow><mrow><mi>q</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mfrac><mrow><msubsup><mi>f</mi><mrow><mi>s</mi><mi>i</mi></mrow><mi>q</mi></msubsup></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>∈</mo><msub><mi>S</mi><mi>q</mi></msub></mrow></munder><mrow><msub><mi>c</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>j</mi><mo>∈</mo><msub><mi>S</mi><mi>q</mi></msub></mrow></munder><mrow><msub><mi>c</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle><mfrac><mrow><msub><mi>p</mi><mi>j</mi></msub></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>i</mi><mo>'</mo></mrow></munder><mrow><msubsup><mi>f</mi><mrow><mi>s</mi><mi>i</mi><mo>'</mo></mrow><mrow><mi>q</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msub><mi>c</mi><mrow><mi>i</mi><mo>'</mo><mi>j</mi></mrow></msub></mrow></mstyle></mrow></mfrac></mrow></math> Where p j is a simulated detected photon number in j -th LOR, f si k is image value of the i -th pixel after the k -th iteration when using projection belonging to s -th subset, and c i,j is an element in the system matrix records the probability that the j -th LOR detects a photon emitted from the i -th pixel. The experiments in this study were based on a practical multi-ring small animal-PET scanner. A Derenzo phantom of F-18 source with 6 groups of hot rods was scanned during experiments. The diameters of the hot rods were 0.75 mm, 1.0 mm, 1.35 mm, 1.7 mm, 2.0 mm and 2.4 mm respectively. In each group, the distance between adjacent rods was twice of the rod diameter. List-mode coincidence event data were collected and histogrammed to a 2D sinogram projection. According to these data we can implement image reconstruction with calculated system matrix for both the polygonal pixel grid and cubic one in this study. Storage size of system matrix, time cost and image quality of can be compared for these two methods. 5 Results and discussion 1 5.1 Structural testing for rotationally symmetric polygonal pixel grid 2 As defined in Eq.(7), three test images are defined with three different segment numbers: P =8, P =9 and P =12. These images are defined with rotationally symmetric polygonal structures, and converted to images expressed with conventional rectangular pixel structure as described in Eq.(5). The pixel size in the converted images is almost 1/5 of that of original images. Three converted images are shown in Fig.11 Fig.11 . No image quality loss is observed after the conversion, which demonstrates the correctness of the image structure conversion process. Figure 10 Figure 10 displays results of the test images defined in Eq.(7) converted to Cartesian coordinate. The segment number P of Figs.10(a-c) is respectively 8, 9 and 12, and the pro-defined images are shown in Fig.6 Fig.6 . The layer-height ρ of the rotationally symmetric structure is 3 mm. The pixel size of the rectangular pixel structure is 1mm, and the final size of the image is 128×128. 5.2 Validation of the generated system matrix through forward-projection process 2 In this test, three different distributions of radioactive sources were defined: 1) A point source located at the center of FOV; 2) A point source located radially 29.6 mm off the center of FOV; 3) A volumetric source which uniformly fills the entire FOV. Three corresponding sinograms were generated by forward projecting the above source images with the generated system matrix, which are shown in Fig.11 Fig.11 . In Fig.11 Fig.11 (a), the sinogram is expected to be a straight line ( r =0), and the generated sinogram agrees well with the expectation. When the point source is 29.6 mm off-centered, the expected sinogram is a sinusoidal curve with 29.6 mm amplitude. In the generated sinogram, the measured amplitude is 29.6 mm. In Fig.11 Fig.11 (c), the covered range in the r direction of the generated sinogram is 29.5 mm. Since the size of FOV is 30 mm, the covered range well matches the FOV size. Therefore, the sinograms generated by system matrix are corresponding to the theoretical expectations, and it approves that the system matrix is calculated correctly. 5.3 Image reconstruction results 2 A hot rod phantom with 6 sections of hot rods (0.75 mm, 1.0 mm, 1.35 mm, 1.7 mm, 2.0 mm and 2.4 mm in diameter, respectively) and 30 mm cross sectional size was scanned using a practical animal PET system. 2D sinograms (108 angular bins and 70 radial bins) were generated from the list-mode event data. The images reconstructed with OS-EM iterations in polygonal pixel grid and cubic grid are shown in Figs.12(a) and (b), respectively. Figure 12(a) Figure 12(a) shows reconstructed images based on a rotationally symmetric polygonal structure, with 12 segments and 79 layers. The layer height was ρ =0.33 mm. After the reconstruction, the images were converted to conventional rectangular pixel structure, with 0.1 mm×0.1 mm pixel size. We can see typical rotational polygonal pixel shaped artifacts in the lower, middle part of the zoomed image section. The right column shows reconstructed images based on a cubic pixel grid with 142×142 pixels of 0.33 mm×0.33 mm sizes. Both in Figs.12(a) and (b), the 1.35 mm hot rod sections can be distinguished. Line profiles in Fig.13 Fig.13 illustrated that the desired 1.35 mm image resolution can be achieved with the proposed rotationally symmetric polygonal structure. Based on this small animal PET system mentioned previously, the memory size of system matrix can be compressed to 46.1 MB when elements are stored as 32-bit float numbers using this polygonal pixel grid. In comparison, when we use the conventional rectangular pixel grid, the required memory size reaches 581.5 MB. The compression ratio of the system matrix size is 12.6:1. It took 2804 s to run 54 iterations in polygonal pixel grid running in a PC with a 2.5 GHz Inter(R) Core(TM) i5-3210M CPU and 4 G memory. The iterations conducted in cubic grid cost 2440 s. Considering the total pixels number in polygonal gird is 1.88 times more than the regular one, the slight increase reconstruction time is acceptable. 6 Conclusion 1 In this work, we designed a novel rotationally symmetric polygonal pixel structure for PET image reconstruction. No image quality loss was observed after the conversion of image in polygonal pixel structure to a conventional rectangular pixel for image display when the area of rectangular pixels is 1/5 of the polygonal ones. Three different distributions of radioactive sources were forward-projected to sonograms and verified the correction of system matrix. Images were successfully reconstructed based on 3D list mode data sinograms of hot rod phantom from a practical small-animal PET. The image resolution of reconstructed phantom image is 1.35 mm. A compress ratio of 12.6:1 of the system matrix size was achieved on this small animal PET, comparing with the conventional rectangular pixel based one. With acceptable time cost, the image resolution based on polygonal pixel structure is as good as the cubic pixel structure based one. In further studies, we will expand our study to the optimization of division mode of rotationally symmetric structure and seek solution for the unevenness of pixel shapes and areas. We can also apply this structure to the 3D system matrix calculation to establish foundations for practical applications in the full 3D PET. References Alessio A , Kinahan P . J Nucl Med , 2006 , 1 : 1 - 22 . J Nucl Med Zhou J , Qi J . Phys Med Biol , 2011 , 56 : 6739 - 6757 . Phys Med Biol Herman G T , Censor Y , Gordon D , et al . J Am Stat Assoc , 1985 , 80 : 22 - 25 . J Am Stat Assoc Kaufman L . IEEE T Med Imaging , 1987 , 6 : 37 - 51 . IEEE T Med Imaging Ansorge R E , Carpenter T A , Clark J , et al . NSS'08 . IEEE , 2008 , 4112 - 4114 . NSS'08 Leroux J D , Thibaudeau C , Lecomte R , et al . NSS'07 . IEEE , 2007, 5 : 3644 - 3648 . NSS'07 Scheins J , Axer M , Pietrzyk U , et al . NSS'07 . IEEE , 2007, 5 : 3458 - 3460 . NSS'07 Scheins J J , Herzog H . NSS'08 . IEEE , 2008 , 4147 - 4150 . NSS'08 Scheins J J , Herzog H . 9th-Fully 3D , 2008 , 49 - 51 . 9th-Fully 3D Scheins J J , Boschen F , Herzog H . IEEE T Med Imaging , 2006 , 25 : 1363 - 1369 . IEEE T Med Imaging Scheins J J , Herzog H , Shah N J . IEEE T Med Imaging , 2011 , 30 : 879 - 892 . IEEE T Med Imaging Cabello J , Rafecas M . Phys Med Biol , 2012 , 57 : 1759 - 1777 . Phys Med Biol Cabello J , Gillam J E , Rafecas M . J Biomed Imaging , 2012 , 5 : 1 - 11 . J Biomed Imaging Siddon R L . Med Phys , 1985 , 12 : 252 - 255 . Med Phys Mora C , Rafecas M . Polar NSS'07. IEEE , 2006 , 5 : 2812 - 2817 . Polar NSS'07. IEEE Shepp L A , Vardi Y . IEEE T Med Imaging , 1982 , 1 : 113 - 122 . IEEE T Med Imaging Lange K , Carson R . J Comput Assist Tomogr , 1984 , 8 : 306 - 316 . J Comput Assist Tomogr Cabello J , Oliver J F , Torres-Espallardo I , et al . NSS'10 . IEEE , 2010 , 3256 - 3261 . NSS'10 10.13538/j.1001-8042/nst.2013.03.011 Supported by Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry and National Natural Science Foundation of China (No.10975086) 2013-03-19 2013-06-01 2021-04-11 2013 Nuclear Science and Techniques 3 24