Athabasca Recon Manual

The following figure shows the defined geometry and coordinate systems, as viewed from above (z and v are up in this figure).

Vector values, written as number tuples, are denoted in the order (x,y,z) for volume, and (u,v) for projections.

Pixels (and their 3D counterparts Voxels) have a finite size, which inevitably leads to the question, "Is the position of pixels denoted by their lower-left corners, or by their centers?" Here we using the center-of-the-pixel convention. Thus the origin of the projections will give the center position of the pixel with index (0,0), and the origin of the reconstruction volume will be the center of the voxel with index(0,0,0). This convention is perhaps more commonly used than the other one when dealing with images. Be aware however, that it has the disadvantage that the numerical value of the origin changes if the pixel/voxel size is changed (as happens if for example you decimate or bin the data to reduce the data size and the resolution). Unfortunately, ImageJ uses a pixel-corner-as-the-origin convention. You’ll just have to add/subtract half a pixel when using ImageJ to compute values such as the center of mass.

The only units that are used in Athabasca Recon are units of length. No particular units are assumed; rather length units are used consistently. Hence if you specify a numerical value for the projection pixel size in microns for example, you must also use microns for the voxel size, and the resulting reconstruction attenuation densities will have units of inverse microns. This may be inconvenient, as for example, the usual units for attenuation density is inverse cm. However, the output can be scaled; see Reconstruction.ScalingFactor.

The execution of athabasca_recon is controlled with a configuration file, which must be specified when it is run, like this:

The configuration file has a simple structure and can be created and edited with any text editor. It can also be given any extension. I prefer to use .conf, but this is not required. Here is a simple example of a configuration file:

The file is divided into sections, denoted with section headings enclosed in square brackets. Within each section is any number of kay/value pairs. The keys and values are case sensitive.

The parameters which can be specified in this file are described in the following sections.In this manual, we will specify the complete key by prepending the section with a period, for example Volume.Dimensions .

For a complete list of possible configuration parameters, refer to Configuration file reference .

Some comments of the formatting follow.

Some values can be specified as a list of multiple values (a "tuple"). There are several acceptable ways to write these. The following are all acceptable and equivalent:

Tuples are ordered as (x,y,z), or (u,v).

Quoting of text values is unsupported, so the following is WRONG:

The following is OK:

Obviously, without quoting and escaping, not every possible string value can be represented. I really don’t expect this to be an insurmountable problem for anyone.

Integer values cannot have a decimal point. The following will cause an error because NumberOfProjections is supposed to be an integer:

Here I’m going to show a simple example of running Athabasca recon.

The input data is the example file projections.mhd, which is distributed with Athabasca Recon. This data file consists of 129 projections of 4 spheres of constant density.

Here is the configuration file we will use. It’s similar to the simple example above, except that we are also writing out the Attenuation Projections, because I want to inspect them. Also, we’re not going to specify anything about the volume; Athabasca Recon chooses pretty reasonable default values.

Before running the reconstruction, we can get the complete configuration, which includes all the default values, as well as values deduced from the data files. The command to obtain the complete configuration is:

The output looks like this:

If you’re satisfied with all the parameters for the reconstruction (and I really do recommend that you make it a habit to look them over before performing the reconstruction), then the reconstruction itself is run with this command:

This will start again with reporting the complete configuration as above. The output that follows is shown below. Comments follow.

After running it, you may want to compare the raw projections with the Attenuation Projections.

The following figure shows the first raw projection, as viewed in ImageJ, scaled to 400%.

This looks like an X-Ray shadow, as it should. Here is the corresponding Attenuation Projection. (If you don’t know what an Attenuation Projection is, see The meaning of Attenuation Projections.)

Notice that the objects in the Attenuation Projection are bright, and the background is dark. The background ought to be zero. We can check that. One way is to use an area selector tool in ImageJ, select some background, and choose Analyze → Measure from the menu. Another visual way is to select the line tool, mark a line across the image, and choose Analyze → Plot Profile. For this data, a plot along the vertical center looks like this:

We see that the background goes nicely to zero.

Now let’s open the reconstructed volume, and scroll to the middle slice (32/64).

We can measure the degree of unevenness in the reconstruction by selecting an area inside one of the spheres with the Circle selection tool (in the image above, the area I selected is the yellow circle). Then from the Analyze menu, select Set Measurements and make sure that Standard Deviation is checked. Finally, select Analyze → Measure. The following results are reported.

We can investigate the sharpness by using the line tool to mark a line passing through the centers of the two larger spheres, and selecting Analyze → Plot Profile. This results in the following plot. The degree of blurring of the sharp edge of the sphere is clearly visible; higher resolution projections would be needed for better sharpness.

As the program processes each batch of projections, it reports its progress through the calculation stages with a series of letters. These letters are flushed to standard out, so they hopefully will appear in your terminal window at the same time as the program sends them. The purpose of these indicators is to give some indication of what fraction of time the program is spending on each calculation step.

It often happens that the center of the projections is not known exactly. Perhaps additionally there is a small unknown misalignment (rotation) between the detector and the axis of rotation. In these cases the Align Projections tool can be used. It is implemented as an ImageJ plug-in.

The Align Projections tool works only on parallel projection data (for now), and it requires that the last projection obtained is at 180º from the first projection (see Projections.ProjectionAt180). It works by aligning the first projection with a mirrored copy of the last projection.

This section will provide a brief tutorial demonstrating the use of the ImageJ plug-in.

The Align Projections tool should be run on Attenuation Projections, and not raw projections (see The meaning of attenuation projections). The reasons for this are:

The following configuration file generates attenuation projections:

Run Athabasca Recon like this to generate the Attenuation Projections:

Open the resulting file attenuation_projections.mhd in ImageJ. Note that, if installed correctly, the MetaImage reader will be found under Plugins → 3D IO → MetaImage Reader…​

You will find the Align Projections tool under Plugins → Align Projections.

Note that the default center pixel is the middle of the projection.

Two things happen when we click the Update button:

As you see from the figure, the red image (which is the first projection) is to the right of the cyan image; hence the actual rotation center is somewhat greater than the trial value for Center Pixel that we used.

Let’s try a larger value for Center Pixel. Try entering 1750 and hitting Update. While we’re at it, let’s reduce the image to the actual region of interest for the reconstruction; enter 200 for Top Border, Bottom Border, and Horizontal Border.

You can see from the image below that the two images are better aligned; the cross-correlation has increased to 0.9928. You may have to repeat this a few times to get close to alignment.

Once you are close to alignment, you can click Optimize and the plug-in will optimize the parmeters by maximizing the cross-correlation. This can take several minutes.

The image below shows the final projection overlay; it appears sharp, with no hint of color, and the cross-correlation has increased to 0.9993 .

Click the Apply to stack and save button. You will be prompted for a file name with extension mhd. In this example we’ll use aligned_projections.mhd. The Align Projections plug-in will now write the entire projection stack to the specified file, with the alignment transformations applied to each image using bi-linear interpolation. This can take several minutes; when finished, the image stack will be opened as a new virtual stack.

To now carry out back-projection, we create a configuration file like the following example. Note that Projections.CenterPixelU does not have to be specified, since the projections in aligned_projections.mhd are centered, which is assumed by default.

The first step in CT reconstruction is converting the raw projections to attenuation values. For the end user, it is useful to be familiar with these attenuation projections, because they can be both inputs and outputs to Athabasca Recon, and they are needed for many refinements (e.g. for correcting for changing beam power). Also, it can be very illustrative to view the Attenuation Projections with an image viewer, particularly if there are problems with the reconstruction.

A little bit of math is useful here, even though we’re not going to go into the whole theory of reconstruction.

A (non-diverging or parallel) ray that passes through a medium decreases in intensity according to

where α is the attenuation coefficient. Integrating this, we obtain the intensity, which is more or less what we measure experimentally:

This is not very practical for direct fast reconstruction, as it is nonlinear in α. However if we take the logarithm, thereby converting it to attenuation, as shown, then we have a nice linear equation relating a quantity, A, derived from the data, to something we want to know, α, the attenuation coefficient, often referred to simply as the density in CT reconstruction, and everything is now copacetic:

To obtain the Attenuation Projections as output from Athabasca Recon, set Output.AttenuationProjectionsFile. To use a set of Attenuation Projections as input, use Input.AttenuationProjectionsFile.

CT images are typically obtained under conditions of non-uniform beam intensity, and real X-ray detectors have dark current, and often a non-uniform response. To deal with these realities, in practice a "dark field" is measured with no X-ray illumination, and then a "bright" field is measured with the X-ray source on at operating level, but no object in the field of view. Then we calculate the attenuation, on a per-pixel basis, according to

In this manual, and in Athabasca Recon, the difference of the bright and the dark fields is referred to as the "flat field".

The files containing the Bright and Dark Fields are specified with Input.DarkFieldFile and Input.BrightFieldFile. The Bright field is required; the Dark field is optional.

Note that noise in the bright and dark fields is additive, as it becomes a systematic perturbation applied to every projection, while noise in an individual projection is effectively decreased by averaging over all projections (by a factor of 1/sqrt(N) ). For this reason, it is desirable to have lower noise in bright and dark fields than in the projections. It is therefore not uncommon to obtain several field measurements so that they can be averaged together. Athabasca Recon will average together multiple bright/dark fields if they are provided.

Before being back-projected, the Attenuation Projections must be convolved with a so-called ramp function. This is typically done by processing with a FFT (Fast Fourier Transform), hence the name of the method "filtered backprojection".

As with the Attenuation Projections, the Filtered Projections can be specified as both input and output; see Input.FilteredProjectionsFile and Output.FilteredProjectionsFile. Unlike the attenuation projections, there are few calculations that can be performed on the filtered projections, and it is rarely illuminating to view them with an image viewer. However, saving the filtered projections can in certain cases speed up processing, as they can, like the Attenuation Projections, be specified as input.

X-Ray detectors sometimes have a few bad pixels. Although a small number of pixels in a high resolution detector represent a negligible loss of information, we do still have to do something explicit with them. In the best base, this avoids ring artefacts, which can arise when bad pixel values are far out of the range of neighbouring pixels. In the worst case, it avoids catastrophic failure of the calculation which can occur when illegal values (such as negative inputs to the log function) result in NaNs which the subsequent FFT then propagates to the entire projection row.

In practice, the actual scheme for dealing with bad pixels matters little as the number of bad pixels is very small and the information content correspondingly so. (If the number of bad pixels is not very small, no miraculous transformation will fix the data and replacement of the detector is required.) Generally, we want to take some sort of average of neighbouring good pixels. The slight complication is that we cannot assume a given neighbouring pixel is good; in fact is not at all atypical in some detectors for bad pixels to occur in clumps. But this is merely a book-keeping issue for the software.

Bad pixels can typically be identified as pixels with excessively large values in the dark field, or pixels with excessively small values in the flat field.

Bad Pixel Correction can be turned on with the setting Reconstruction.BadPixelCorrection. This turns on both automatic identification of bad pixels, as well as automatic correction/fudging. See also Reconstruction.FlatFieldBadThreshold and Reconstruction.DarkFieldBadThreshold.

For synchrotron tomography, the beam power decays with time, and this can be significant over the time of the measurement.

There are several possible approaches to correcting for this, which are described below. The method is selected with Reconstruction.BeamPowerCorrection

CT Reconstruction has the characteristic of amplifying high frequency noise. Without additional processing, reconstructed volumes typically have a speckled appearance, with a large amount of high-frequency noise. Of course, this can be dealt with by post-processing with additional software. However, it is common to employ a low-pass filter in the reconstruction process.

The options for low pass filtering are described below. The filter is selected with the parameter Reconstruction.SmoothingFilter.

Back-projection is essentially ray-tracing from each voxel to the projection. Of course, the rays don’t necessarily hit exactly in the middle of a projection pixel, so some kind of interpolation scheme is required. There are a couple of possible approaches. These are selected with the parameter Reconstruction.PixelInterpolation. The options are discussed below.

The value of the nearest pixel is used. This is the fastest and least accurate method.

A bi-linear interpolation of the 4 nearest pixel centers is used. This is slower than nearest neighbor, but more accurate.

More theoretically accurate schemes exist, such as bi-cubic, but they generally have the property of amplifying noise, and are not supported by Athabasca Recon.

One drawback of bilinear interpolation is that the domain over which interpolation can be performed is limited to the interior extents of the projection; that is, to the rectangular region with corners at the pixel centers of the corner pixels. In contrast, nearest-neighbor interpolation is valid over a domain going right to the edges of the outer pixels (the exterior extents). Thus the domain of nearest-neighbor interpolation is 1/2 pixel wider on every side. In practice, this can make a difference if you want to reconstruct a volume that goes right to the edge of a projection row. It sometimes happens that, using bilinear interpolation, the outer slices coincident with the top and bottom projection rows are calculated as identically zero. To avoid this problem, use the method bi-linear with fallback, which uses bilinear interpolation inside the interior extents, and nearest neighbor interpolation within the 1/2 pixel wide band outside of the interior extents.

Some care is required when doing reconstructions with large voxels (i.e. low volume resolution). By "large" I mean as compared with the pixel spacing.

To see why, consider a voxel that has an edge length 4 times as big as the projection pixel spacing. Conceptually, each voxel, when back-projected onto the projection, would cast a shadow over roughly 4 × 4 = 16 pixels (actually more, since the voxel can be rotated at some angle). But the pixel interpolation methods currently available in this program are nearest neighbor on the central ray (uses data from 1 pixel) and bilinear interpolation on the central ray (uses data from 4 pixels, although strictly its information content is only 2 pixels). Thus, assuming bilinear interpolation, only about 1/8 of the available information is used in the reconstruction. This has the usual undesirable implications for noise and accuracy.

An ideal reconstruction program, if asked to reconstruct a volume with a voxel size much larger than the pixel spacing, would first down-sample the projection data, using a summing or averaging procedure. Athabasca Recon is not currently that program.

The work-around, if you want high-quality low-resolution reconstructions, is to pre-downsample the projections using some other program.

In this section the input files are specified. You need to specify exactly one of RawProjectionsDataFile, AttenuationProjectionsFile and FilteredProjectionsFile. (Note that there are similar parameters in the output section, which are independent of the input parameters.)

The input data file of unprocessed projections.

The type of file is automatically determined by the extension.

The data file containing the Dark Field.

The data can consist of a single dark field (i.e. 2D data), or any number of dark field measurements (i.e. 3D data). In the latter case all the available fields are averaged together with floating point precision.

Some CT file formats store this integrally, in which case it is not necessary to specify a value for this parameter.

The Dark Field is optional; it can be left unspecified if no Dark Field is available.

The data file containing the Bright Field.

The data can consist of a single bright field (i.e. 2D data), or any number of bright field measurements (i.e. 3D data). In the latter case all the available fields are averaged together with floating point precision.

Some CT file formats store this integrally, in which case it is not necessary to specify a value for this parameter.

A Bright Field is required.

The data file containing the bright field as measured after the scan. This data is only required if using the BeforeAndAfterBrightField setting for Reconstruction.BeamPowerCorrection.

The input data file of projections, converted previously to attenuation values. This setting is typically useful if you want to perform only the attenuation calculation at one time (saving the result with Output.AttenuationProjectionsFile), and then re-use those when running the program later.

The type of file is automatically determined by the extension.

The input data file of projections, converted previously to attenutation values and DFT filtered. This setting is typically useful if you want to prepare the projections for back-projection (saving the result with Output.FilteredProjectionsFile), but not actually perform the back-projection at the same time. The back-projection alone can then be done using this input file parameter.

The type of file is automatically determined by the extension.

The name of the file to write the reconstructed volume to. At the moment, only .mhd output files are supported.

If the parameter is not present, then no back-projection will be performed, although projection processing and filtering will proceed as usual (and will be saved to the specified files if you choose).

If specified, then the attenuation projections will written to the specified file. At the moment, only .mhd output files are supported.

If specified, then the FFT filtered projections will written to the specified file. At the moment, only .mhd output files are supported.

If specified, and Reconstruction.BeamPowerCorrection is not None, then the attenuation correction for each projection will be written to the specified file (in text format, suitable for import into a spreadsheet).

In addition, a best linear fit to the logarithmic beam power corrections is automatically performed and reported.

These are parameters which describe the projection data. For many types of projection data files, several of these parameters can be read from the data file itself. In such cases, it is not necessary to also specify them in the configuration file. If however you do so, the value in the configuration file overrides the values read from the projection data file.

Specifies the data type of the projection data. Valid values are INT8, UINT8, INT16, UINT16, INT32, UINT32, FLOAT32, FLOAT64 .

Specifies the number of pixel dimensions of the projections. Must be a length 2 tuple in the form (dim u, dim v).

Specifies the number of projections.

Specifies whether the number of projections includes the last projection at 180º or not. Valid values are True and False. If False, the angle increment is 180º/NumberOfProjections and the last projection is back-projected. If True, the angle increment is 180º/(NumberOfProjections + 1) and the last projection is not back-projected.

The default is True.

Specifies the spacing of the pixels on the projections. Units are real-space units (e.g. mm, or whatever you prefer to have the output in). Must be a length 2 tuple in the form (pixel spacing in the u direction, pixel spacing in the v direction).

If this parameter is not specified, and the pixel dimensions cannot be determined from the input files, then (1,1) will be assumed.

The projection of the center of rotation on the projections. Units are pixels, as an offset from the u origin of the projection.

The default is the center of the projection (i.e. (N_U-1)/2, where N_U is the number of pixels in the U direction).

Specifies the offset of the projections in the v direction in real units; in other words the location in the z direction of row 0 (or of v=0).

The default is to center the projections around z=0.

Specifies that the projections were obtained in a direction opposite to the standard convention for this software. See Geometry and Coordinate System .

This parameter can be used to process only every n^th projection. For example, if set to 2, only every second projection in the input file will be used.

The volume parameters describe the volume which is to be reconstructed. These parameters are superfluous if Output.VolumeFile is not set.

The dimensions of the reconstruction volume in voxels, specified as a tuple in the form (dim x, dim y, dim z).

The spacing of the reconstruction volume voxels, specified as a tuple in the form (spacing x, spacing y, spacing z).

The default is the same size as the projection pixels, provided that the pixels are square. If the pixels are non-square, there is no default voxel size.

The origin, in real space, of the reconstruction volume (i.e. the location in real space of the center of voxel 0,0,0). Should be specified as a tuple in the form (origin x, origin y, origin z).

If not specified, the volume will be centred at the origin of the coordinate system.

These parameters modify how the reconstruction is performed.

Sets a scaling factor for the reconstructed volume. The default is 1.

As an example, if you specify the projection pixel spacing and the voxel size in microns, you likely want a ScalingFactor of 1E4 in order to obtain attenuation densities in inverse cm.

Sets the method for correcting bad pixels. Valid values are None and Averaging. The Averaging correction will take the average value of the 4 nearest neighbors, or, if any of the nearest neighbors are themselves marked as bad pixels, will take the 4 nearest good pixels to determine an average. See the discussion section Bad pixel correction.

If Bad Pixel Correction is being used, sets the threshold below which pixels in the flat field are flagged as bad. The default is 10. Note that this is an integral value which is reasonably plausible for discrete-valued input data, as raw projection data typically is. (i.e. If your detector count is below 10 in the flat field, something is wrong.) However, it is likely to be totally non-sensical for scaled floating-point input data.

If Bad Pixel Correction is being used, sets the threshold above which pixels in the dark field are flagged as bad. The default is for this value not to be set, in which case automatic detection of bad pixels from the dark field is disabled.

Sets the method for correcting for beam power decay. See Correction for decaying beam power. Valid values are shown in the following table.

Set the independent variable for the beam power correction. Valid values are ProjectionNumber and Time. The default is ProjectionNumber.

Note that ProjectionNumber is zero-indexed; the first projection is projection number 0.

This value sets the attenuation correction. The constant term is subtracted from each attenuation projection.

Note that this value is typically calculated as the logarithm of a ratio of beam powers or intensities.

See Correction for decaying beam power for details.

This value sets the attenuation correction. The linear term is subtracted from each attenuation projection after being multiplied by the projection index (which is 0 for the first projection), or by the time, depending on the value of BeamPowerIndependentVariable. The manual method uses both these values; some other methods allow only the constant term to be independently specified.

See Correction for decaying beam power for details.

This setting is used only for BeamPowerCorrection = NullProjectionEdge. It determines the width (in pixels) of the strip on either side of the projection that is used to determine the background beam power. This strip should never have any object causing attenuation in it, and should be illuminated by the beam.

Sets the method for interpolating the pixel values of the projections. Valid values are NearestNeighbor, Bilinear and BilinearWithFallback. The default is BilinearWithFallback. NearestNeighbor will give a speed increase, at the cost of some accuracy.

See Pixel interpolation.

Selects a smoothing filter to apply to the projections before back projecting. See Low pass filtering. The default is Gaussian.

Applies only when Reconstruction.SmoothingFilter = Gaussian . Sets the σ value of the convolution gaussian, in real space, in units of pixels. The default is 0.5 .

Sets the smoothing filter frequencies parameters. Only relevant for appropriate filter types. See Low-pass filtering. No default; must be set if required.

These parameters affect how the software runs. They should not affect the results, except indirectly by for example changing order at which data is processed (which can affect round off errors).

Current options are SingleThreaded and MultiThreaded. The default is MultiThreaded. As SingleThreaded is slower but produces the same results, it is primarily intended to be used for debugging, since it can be difficult to debug a multi-threaded program.

This option only applies when Engine is MultiThreaded. Valid values are any positive integer, and Automatic. Automatic will select a value equal to the number of CPU cores in the system. Automatic is the default.

Specifies the maximum memory that will be used to store volume data. This will determine how many passes are required to reconstruct the entire volume. Larger values (fewer passes) are typically faster, however if set too large, swapping of virtual memory to disk will occur, which will result in very slow calculation times. Valid values are Automatic (which is the default) or a numerical value appended with MB or GB. For Automatic, the value used will be the installed system physical minus 1GB.

Filtering module selects the module to use for convolution/filtering. The choices are as shown in the following table.

ImageJ can be obtained from http://rsbweb.nih.gov/ij/download.html . Version 1.45 or newer is required.

The Align Projections plug-in is distributed with Athabasca Recon. Simply put the Align_Projections.jar file in the plugins directory of the ImageJ.

A few other plug-ins are required or recommended.

Align Projections requires Apache Commons Math, which can be obtained from http://commons.apache.org/math/download_math.cgi . Again, just put the commons-math-2.2.jar file in the plugins directory of ImageJ.

The plug-in for reading and writing ITK MetaImage files is available at http://ij-plugins.sourceforge.net/ . The one you want is "ij-Plugins Toolkit".

The Plugin for reading Hamamatsu’s SimplePCI files (.cxd) is available from http://www.loci.wisc.edu/bio-formats/imagej .

In the current version Athabasca recon only reads ITK MetaImage files (and raw data files). This format was chosen because:

ITK MetaImage files actually consist of two files; one is a raw data file (extension .raw) and one is a simple text file that describes the data (extension .mhd). The text file can be created or modified with any text editor.

Here is an example .mhd file. Given a raw file with known parameters, you can construct the corresponding .mhd file using this template.

To open the .cxd file, you have to to go the "Plugin" menu as shown, not to the "File" menu.

The important option is "Use virtual stack". This allows you to open very large data files without requiring a lot of RAM. You should also select "View stack with: Standard ImageJ".

When writing out files as MetaImage, make sure that the option "Save in single file" is unselected.

A compiler. Athabasca Recon has been compiled with gcc on Linux, XCode 6.4 on OS X, and with Visual Studio 2013 on Windows. Other compilers are also likely to work.

CMake. CMake is used to support cross-platform building of Athabasca Recon. You can get cmake from http://www.cmake.org/ .

boost. Boost is a high quality collection of portable C++ libraries. You can get boost from http://www.boost.org/ .

An FFT library. On OS X, Apple’s vDSP is used by default; it is already present on any OS X system, and requires no additional installation. The other FFT library that Athabasca Recon supports is FFTW. On Linux systems, you will typically use your distro’s package management tool to install FFTW. For Windows users, see http://www.fftw.org/install/windows.html . Make sure that you follow the instructions there about running lib with the /machine:x64 option. In any case see the example builds below.

Google Test. Google Test is used for the unit tests. It is possible to build Athabasca Recon without Google Test by setting ENABLE_TESTING to OFF. However if you are modifying the source code, I really recommend that you build and run the unit tests. You can get Google Test from https://github.com/google/googletest .

1. Get FFTW (Linux only)

If you are on OS X, skip this step, since Apple’s vDSP will be used. How you get fftw will depend on the Linux distro you are using.

On RedHat Enterprise Linux or CentOS, the following should work:

On Ubuntu Linux, the following should work:

2. Get and compile boost.

Here is how I compile boost.

3. Get and compile Google Test.

Unpack the source code. On Linux and OS X the defaults are all good, so just run

4. Create a build directory.

From the source code directory of Athabasca Recon, do

5. Run cmake.

Before running CMake, we are going to set a environment variable to help CMake find boost. This just simplifies things a bit.

Of course, this is just an example: you have to specify the actual path to boost.

Run ccmake from the build directory, specifying the source directory on the command line:

Hit c for configure.

6. Set CMake variables.

Set CMAKE_BUILD_TYPE to Release (or Debug, if you want to be able to debug the program; but it runs slower).

We need to specify the location of Google Test, so hit t for advanced settings, and modify the following values. Of course you have to use paths appropriate for your system.

Currently, in Linux we also have to add -fpermissive to CMAKE_CXX_FLAGS. This is an advanced setting, which you get to by hitting t.

Hit c again for configure. Now hit g for generate. CMake will exit.

7. Run make.

That’s it. You should now have an executable athabasca_recon.

1. Get FFTW

As mentioned above, download the DLLs from http://www.fftw.org/install/windows.html . Unpack the zip files anywhere. For this example, I used C:\Users\Eric\Install\fftw-3.3.4 .

Now, open a VS2013 Native Tools Command Prompt, change to the folder where you unpacked fftw and run

2. Get and compile boost.

Here is how I compile boost.

3. Build Google Test (optional)

https://github.com/google/googletest . Instructions are at https://github.com/google/googletest/blob/master/googletest/README.md , but they build be default 32 bit libraries, so we have to modify a bit.

Unpack the zip file somewhere. Start a VS2013 x64 Native Tools Command Prompt. cd to the directory you unpacked the gtest source code, then

Hit Configure. You’ll get a pop-up window; choose your development environment ("Visual Studio 12 2013 Win64" is what we use).

It is easiest to link to Google Test dynamic libraries, so select BUILD_SHARED_LIBS. Hit Configure, then Generate. Close CMake. In the build directory is a Visual Studio solution file that you can now open to build Google Test. Make sure to build the Release configuration.

4. Start CMake.

It is best to launch CMake from a VS2013 Native Tools Command Prompt.

First change to the directory where you have the Athabasca Recon source code:

or wherever it is. Now create a build directory

And we will set an environment variable to help it find boost.

Of course you have to specify the folder where boost is on your system. Yes, those are forward slashes. Backwards slashes work in principle as well.

Now we are ready to launch CMake:

5. Hit Configure.

You’ll get a pop-up window; choose your development environment ("Visual Studio 12 2013 Win64" is what we use).

Click Finish. You’ll get an error message: "Error in configuration process, project files may be invalid". This is normal.

6. Identify the location of Google Test and FFTW

If building with Google Test, then specify the location of the libraries and the include directory, as shown in the screenshot. Otherwise, set ENABLE_TESTING to OFF.

For FFTW you need to manually specify both the include path and the two required libraries, as shown in the screenshot.

Now hit Configure again. This time there should be no error messages, and you should get the message "Configuring done".

7. Hit Generate.

You will now have a Visual Studio solution file, Athabasca_Recon.sln, in your build directory. Double click on it to open it in Visual Studio.

8. Build.

Select Release from the Solution Configurations drop-down box. Select the project ALL_BUILD in the solution explorer, then from the menu select Build → Build Solution.

9. Make sure that the FFTW libraries can be found.

You will need to make sure that either the path to the FFTW libraries is added to your PATH variable, or you can simply copy the FFTW .dll files to the build\Release directory.

If you’re modified the source code, it is important to run the tests to verify the code. Unit test coverage is not complete, but it catches many errors.

Before running the tests, make sure that all necessary dynamically linked libraries can be found.

For example, on Windows:

On OS X:

On Linux it is gnerally not necessary to set any library search paths, as we link google test statically and FFTW is installed in system locations by the package manager. (However should it be necessary, we could add a path to LD_LIBRARY_PATH.)

The tests can be run from the build directory by running

The output will look something like this:

Each of the tests listed is actually a collection of Google Test units tests. If a test fails, you can examine the detailed output in the file Testing/Temporary/LastTest.log (in the build directory).

You can also run each Google Test suite individually and with more verbose output, for example

Theoretical derivation of filtered back-projection shows that convolution with a ramp function is required, which in continuous k-space is

Convolution is usually performed as multiplication in k-space using FFTs, as this is much the fastest way to perform this operation. There are however, two potential pitfalls. The first is that the above function fails to exist in continuous real space. That at first might appear not to affect us so much as we are doing a discrete transform anyway. The discrete real space filter function is

Now of course we can only deal with discrete functions with finite support. This seems a bit obscure and theoretical, but the practical implication is that Discrete Fourier Transform of g_i defined for i < N is not longer exactly |k|. The correct procedure is to perform a DFT of the given g_i to obtain the correct G. [1]

The other thing to be careful of is that when using an FFT to perform a convolution, we must zero-pad (in real-space) the input rows out to at least twice the original length in order to avoid wrap-around, which in this application would be an error. Note that that ramp function itself should not be zero-padded, but you must use the length equal to the padded length.

One more thought: It may have occurred to you to wonder whether you can get away with a real-space kernel for the ramp function that is shorter than the projection row length and hence potentially a shorter FFT. (Well, it occurred to me.) The answer is no. Basically because the real-space kernel has a 1/r dependence, and so its absolute sum diverges; given a desired accuracy, there is nowhere safe to truncate it. A couple of consequences are: