Reconstruct an image from C++

In this tutorial you will write a short C++ program that uses the FORGE library directly. It builds an encoding operator, simulates k-space data from a known image, and reconstructs the image with FORGE's conjugate gradient solver.

You need a FORGE source checkout that builds, from Build and install FORGE. The library headers are not installed, so in this tutorial you build your program inside the FORGE build.

Write the program

In the top of your FORGE checkout, make a directory for the program:

mkdir myrecon

Save this as myrecon/main.cpp:

#include "Core/forge.h"
#include "Operators/Gdft.h"
#include "Penalties/QuadPenalty.h"
#include "Solvers/solve_pwls_pcg.hpp"

#include <iostream>

using namespace arma;

int main()
{
    const uword Nx = 16, Ny = 16;
    const uword nImage = Nx * Ny; // image pixels
    const uword nData = Nx * Ny;  // k-space samples, one per pixel

    // Pixel positions (fraction of the field of view) and a Cartesian
    // trajectory (cycles per field of view). y varies fastest, as in
    // every FORGE image.
    Col<float> ix(nImage), iy(nImage), iz(nImage, fill::zeros);
    Col<float> kx(nData), ky(nData), kz(nData, fill::zeros);
    for (uword x = 0; x < Nx; x++) {
        for (uword y = 0; y < Ny; y++) {
            uword n = y + x * Ny;
            ix(n) = float(x) / Nx;
            iy(n) = float(y) / Ny;
            kx(n) = float(x) - Nx / 2.0f;
            ky(n) = float(y) - Ny / 2.0f;
        }
    }

    // No off-resonance, so the field map and sample times are zero.
    Col<float> fieldMap(nImage, fill::zeros); // rad/s
    Col<float> times(nData, fill::zeros);     // s

    Gdft<float> G(nData, nImage, kx, ky, kz, ix, iy, iz, fieldMap, times);

    // A square in the middle of the image.
    Col<cx_float> truth(nImage, fill::zeros);
    for (uword x = 4; x < 12; x++)
        for (uword y = 4; y < 12; y++)
            truth(y + x * Ny) = cx_float(1.0f, 0.0f);

    // Simulate the scan with the forward transform.
    Col<cx_float> kspace = G * truth;

    // Reconstruct with 20 conjugate gradient iterations.
    forgeCol<forgeComplex<float>> y(kspace);
    forgeCol<float> W(Col<float>(nData, fill::ones));
    forgeCol<forgeComplex<float>> x0(Col<cx_float>(nImage, fill::zeros));
    QuadPenalty<float> R(Nx, Ny, 1, 0.0f);

    forgeCol<forgeComplex<float>> xhat = solve_pwls_pcg<float>(x0, G, W, y, R, 20);

    float err = norm(xhat.getArma() - truth) / norm(truth);
    std::cout << "relative error: " << err << std::endl;
    return 0;
}

The program does four things:

1. It describes a 16 × 16 image and a fully sampled Cartesian trajectory. Pixel positions are fractions of the field of view, and k-space positions are in cycles per field of view, from −N/2 to N/2.
2. It builds a Gdft, the exact non-uniform DFT. The field map and sample times are zero, so there is no off-resonance.
3. It makes k-space data with the forward transform, G * truth.
4. It hands the data to solve_pwls_pcg and compares the result with the true image.

Build it

Add two lines to the end of FORGE's top-level CMakeLists.txt:

add_executable(myrecon myrecon/main.cpp)
target_link_libraries(myrecon PRIVATE ForgeCommon_cpu)

Linking ForgeCommon_cpu brings in FORGE's include paths and its dependencies. Then build:

cmake -B build
cmake --build build --target myrecon

Run it

./build/myrecon 2>/dev/null

You see something like this:

[info] [solve_pwls_pcg.hpp:436] Starting PCG solver: 20 iterations requested
[info] [solve_pwls_pcg.hpp:483] Terminating early: zero gradient at iteration 2
relative error: 2.85696e-07

The exact digits depend on your machine. The solver stopped after 2 iterations instead of 20: a fully sampled Cartesian DFT is so well conditioned that the solver reaches the answer almost at once. Real non-Cartesian, undersampled problems need more iterations, as you saw in Your first reconstruction.

The 2>/dev/null hides the solver's progress reports, which it writes to stderr as JSON lines, one per iteration. JSON lines protocol describes them.

What you did

You built an encoding operator, used its forward transform to simulate a scan, and recovered the image with the solver that the command-line programs use.

From here:

• Operators explains the operator convention and how operators nest, for example SENSE around TimeSegmentation around Gnufft.
• forgeCol and forgeMat describes the vector types the solver works on.
• The C++ API reference lists every class.
• forge/Tests/operatorTests.cpp has more examples, including Gnufft and SENSE.