Next: Samples Up: Walk Through of the Previous: OSEM Style files

Weight Matrix Definition

The PBD matrix or weight matrix is the main complexity in the implementation of an OSEM algorithm. If it is calculated on the fly it slows the iterations down dramatically (more so in the case of attenuation), also it complicates the coding greatly. If it is precalculated in full it is too large to store. Dealing with the various trade offs is an interesting design process. The methods used within the R program (and its counterparts) are detailed below. In the discussion N refers to the size of the reconstruction. e.g. a reconstruction of pixels would have size N=64. The weight matrix is divided into three separate sections, being Geometry, Point Spread Function and Attenuation. Each section forms a lookup table that allows accurate and fast location of weight for each pixels' projections.

Geometry

The Geometry describes for each pixel and projection, the projection bin position that forms a normal from that projection. In the case of the -FanBeam option the projection bin position is adjusted accordingly. Also the distance from the pixel to the camera is calculated. Geometry does not vary from slice to slice, therefore we only need to calculate for a single slice.

There are pixels in one slice and N projections, therefore there are pieces of information in the geometry lookup table. Size of a piece of information is usually 6 bytes. e.g. N=128 yields Geometry = 12Mb.

Point Spread Function

The Point Spread Function describes the spread of the counts about the geometric projection point. The spread is based on the idea of camera blur of a normal distribution. The actual distribution is created using the -CollConst and -CollScale options to produce a lookup table for each distance from the camera and each offset within the projection bin. Pixels' projections have no offset along the slice axis, therefore spread through slices is symmetrical.

Formula for point spread is , where and d is distance from center of geometric projection in mm. PSxWidth and PSzWidth will never cover the entire domain of the ditribution, so to ensure distribution sums to 1.0 the used part of the distribution is scaled accordingly.

Table size is not related to N, but instead to -PSxWidth and PSzWidth options. Also scales according to physical size of reconstruction, usually 400 mm. Size of table is . e.g. By default PSxWidth=1 and PSzWidth=0, so size of table is 9000 elements each of size 4 bytes yielding 36000 bytes.

Attenuation Coefficients

The Attenuation coefficients scale the projected values by number in range 0.0..1.0. Attenuation coefficient is assumed to be along the path normal to the camera from a pixel. To save memory, the values are stored in a single byte and scaled at run time to correct range. If attenuation is not used, table is not created.

Formula for AC is .

There are pixels and N camera positions for which to store a value, therefore the table size is . Each element is a single byte. e.g. N=128 yields Attenuation Coefficients Table = 256Mb. Therefore you may wish to use some approximation. eg. use every 4th slice of the attenuation mu map, this would yield 64Mb. The -Attenuation option lists your methods of dealing with the large memory usage problems of attenuation correction.

Scatter Loss Coefficients

The Scatter coefficients are quite similar to the Attenuation coefficients. They scale the projected values by number in the range 0.0..1.0. Scatter coefficient is assumed to be along the path normal to the camera from a pixel. To save memory, the values are stored in a single byte and scaled at run time to correct range. If scatter is not used, table is not created.

Formula for SC is

There are pixels and N camera positions for which to store a value, therefore the table size is . Each element is a single byte. e.g. N=128 yields Scatter Coefficients Table = 256Mb.

Scatter Blurring Affect

This affect is modelled using fourier techniques in a way separate from the PBD matrix calculation. A kernel (response function) is calculated based upon value of -ScatSd and -ScatSs. , where s is slice distance, d is detector distance.

The scatter projected data is convolved with the kernel, then combined with normal projection data according to the ratio (1-f).N + f.S where f is the value of -ScatFract, N is normal projection, S is scatter projection.

The Look Up

Given a pixel and a projection we look up the Geometry table to find the projection bin and offset within the projection bin that forms the normal to the projection, and the distance of pixel from the camera. Given the distance and offset we lookup the Point Spread Function that gives us a matrix of projection values totaling 1.0. We now scale that matrix by look up from the Attenuation Coefficients table, indexed by pixel and projection. We now scale that matrix by look up from the Scatter Coefficients table, indexed by pixel and projection. We end up with a pixel, and its relevant weights for given projection.

Given final matrix has M elements, cost is M+1 multiplications and less than a dozen table dimensions indexed.