Logo C++ Online Judge

C++

时间限制:3 s 空间限制:512 MB
TA:

2025 Fall Practice 2

A. Computed Tomography 🔴 (Hard)

统计

1. Background

In Computed Tomography (CT), we want to reconstruct the internal structure of a 3D object from many X-ray projections taken at different angles.

In practice, the scanner does not reconstruct the whole 3D volume at once. Instead, it scans the patient slice by slice along the body axis. For each fixed position of the table:

  • The X-ray source and detector rotate around the patient.
  • Many 2D projections of that single cross-sectional slice are acquired at different angles.
  • From these projections, the image of that slice is reconstructed.

Repeating this process at many adjacent positions of the table gives a stack of slices, which together form the 3D volume:

Each frame in the animation below is one such reconstructed slice:

For understanding the reconstruction algorithm, it is enough to focus on one slice at a time. We can idealize the problem as purely 2D: given many 2D projections of a single cross-section, how do we reconstruct the 2D density inside that slice? Details such as table motion or the stepping between slices do not affect this basic 2D reconstruction model.

In a simple 2D model:

  • The object is represented by an $n \times n$ grid (matrix).
  • Each cell $(i, j)$ stores a density value.
  • Each X-ray is modeled as a straight ray passing through this grid.
  • The projection value of one ray is the (approximate) line integral of the density along that ray.

In a discrete grid, this is approximated as:

$$ \text{projection} \approx \sum_{i,j} \bigl(\text{density}_{i,j} \times \text{length of ray inside cell }(i,j)\bigr). $$

The provided C++ program implements a simple iterative reconstruction scheme:

  1. It reads projection data from standard input.
  2. For each ray (given by an angle and an offset), it:
    • Computes the current predicted projection from the current matrix (forward projection).
    • Compares it with the measured projection.
    • Updates the matrix along that ray to reduce the error.

You do not need to know CT reconstruction theory in detail, but you should understand basic 2D geometry and C++.

2. Required Knowledge

This problem focuses on:

  1. Code reading and function usage
    • Understanding a small C++ class with helper functions.
    • Calling existing helper functions with correct arguments at the right places.
  2. Basic C++
    • struct and member functions.
    • std::vector as a 2D matrix.
    • Standard input/output with std::cin and std::cout.
  3. Elementary 2D geometry / graphics
    • Representing 2D points as $(x, y)$.
    • Basic vector operations: addition, scalar multiplication.
    • Parametric form of a line segment: $$ P(t) = \text{start} + t \cdot \text{dir}, \quad 0 \le t \le 1. $$
    • Euclidean length of a vector: $$ \|\text{dir}\| = \sqrt{\text{dir}_x^2 + \text{dir}_y^2}. $$
    • Intersection parameters $t$ with vertical and horizontal grid lines.

3. Role of computeProjection

The function computeProjection performs the forward projection:

  • Input: an angle angle (in radians) and an offset defining one ray.
  • It constructs the ray’s start and end points based on this angle and offset.
  • For each grid cell, it:
    • Uses getRayLength to compute how far the ray travels inside that cell.
    • Multiplies this length by the cell’s density and accumulates the result.

Thus, for the current matrix, computeProjection tells us what projection value we expect along this ray. In the reconstruction loop, this predicted value is compared with the measured projection, and the difference is used to update the matrix.

4. Input and Output

Input (from stdin)

  1. An integer $n$: the matrix size (the image is $n \times n$).
  2. An integer $k$: the number of projection angles.
  3. A line with $k$ floating-point numbers: the projection angles in degrees.
  4. For each angle index $a = 0, 1, \dots, k-1$:
    • An integer $m_a$: number of rays for this angle.
    • Then $m_a$ floating-point numbers: the measured projection values for that angle.

Output

After running a fixed number of iterations (50 in the code), the program prints the reconstructed matrix: one row per line, elements separated by spaces, with two digits after the decimal point.

The output is not unique. Please output any of them that satisfies the projections.

Example

Input:

4
2
0 90
8
0 0 8.0 16.0 5.2 5.4 0 0
8
0 0 5.4 5.2 16.02 7.98 0 0

Output:

1.98   3.96   1.02   1.04
3.94   7.88   2.06   2.12
1.02   2.06   1.04   1.08
1.04   2.12   1.08   1.16

The projections in the input are from two angles, and there are four offsets in each direction:

     |      |      |      |
     |      |      |      |
 -- 1.98   3.96   1.02   1.04 -->  8.00
 -- 3.94   7.88   2.06   2.12 --> 16.00
 -- 1.02   2.06   1.04   1.08 -->  5.20
 -- 1.04   2.12   1.08   1.16 -->  5.40
     |      |      |      |
     v      v      v      v
    7.98  16.02   5.20   5.40

Example 2

Input:

6
5
0.00000000 22.50000000 45.00000000 67.50000000 90.00000000
12
0.0000000000 0.0000000000 0.0000000000 11.7070000000 11.8320000000 10.8890000000 14.9480000000 12.3630000000 13.6360000000 0.0000000000 0.0000000000 0.0000000000
12
0.0000000000 0.0000000000 0.8817005863 5.1218104808 10.1265454193 14.8563922416 16.3029107738 13.5070719121 8.3002717584 0.5920112696 0.0000000000 0.0000000000
12
0.0000000000 0.0000000000 2.6070263730 3.7775127567 10.8919419526 13.8706412447 16.2771866993 12.9005404469 6.4544913125 1.7504683755 0.0000000000 0.0000000000
12
0.0000000000 0.0000000000 0.8817005863 4.1834712680 13.8281520182 12.1469586623 13.9438717536 15.7165407402 6.8324796860 0.5920112696 0.0000000000 0.0000000000
12
0.0000000000 0.0000000000 0.0000000000 12.3720000000 10.4750000000 11.7010000000 14.5750000000 13.1680000000 13.0840000000 0.0000000000 0.0000000000 0.0000000000

Output:

2.013000 1.180000 2.608000 1.311000 1.524000 3.071000
3.144000 2.608000 1.542000 1.710000 1.861000 0.967000
2.682000 2.490000 1.779000 1.052000 2.308000 0.578000
2.185000 1.755000 2.496000 3.403000 2.223000 2.886000
0.998000 1.995000 3.214000 1.835000 1.648000 2.673000
2.062000 3.140000 2.936000 2.390000 0.911000 2.197000

5. Your Task

Fill in all places marked with TODO in the code below.

You are not allowed to change other places in the code.

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

const double PI  = 3.141592653589793;
const double EPS = 1e-10;

struct Point {
    double x, y;
    Point(double x = 0.0, double y = 0.0) : x(x), y(y) {}

    Point operator-(const Point& p) const { return /* TODO */; }
};

// Convert degrees to radians
double deg2rad(double degrees) {
    return /* TODO */;
}

// Euclidean length of a 2D vector
double vectorLength(const Point& v) {
    return /* TODO */;
}

// Intersection parameter t when ray crosses vertical line x = boundaryX
double intersectionParamX(double boundaryX, const Point& start, const Point& dir) {
    return /* TODO */;
}

// Intersection parameter t when ray crosses horizontal line y = boundaryY
double intersectionParamY(double boundaryY, const Point& start, const Point& dir) {
    return /* TODO */;
}

class CTReconstructor {
private:
    int n;
    std::vector<std::vector<double>> matrix;
    std::vector<std::vector<double>> projections;
    std::vector<double> angles;

    bool buildRaySegment(double angle, double offset, Point& start, Point& end) const {
        const double INF = 1e100;

        double dx = std::cos(angle);
        double dy = std::sin(angle);
        double nx = -dy;
        double ny = dx;

        double cx = static_cast<double>(n) / 2.0;
        double cy = static_cast<double>(n) / 2.0;

        double px = cx + offset * nx;
        double py = cy + offset * ny;

        double t_min = -INF;
        double t_max = INF;

        auto update = [&](double coord_start, double dir_comp, double lower, double upper) -> bool {
            if (std::abs(dir_comp) < EPS) {
                return coord_start >= lower - EPS && coord_start <= upper + EPS;
            }
            double t1 = (lower - coord_start) / dir_comp;
            double t2 = (upper - coord_start) / dir_comp;
            if (t1 > t2) {
                std::swap(t1, t2);
            }
            t_min = std::max(t_min, t1);
            t_max = std::min(t_max, t2);
            return t_max >= t_min - EPS;
        };

        if (!update(px, dx, 0.0, static_cast<double>(n))) {
            return false;
        }
        if (!update(py, dy, 0.0, static_cast<double>(n))) {
            return false;
        }
        if (t_max < t_min) {
            return false;
        }

        start = Point(px + t_min * dx, py + t_min * dy);
        end   = Point(px + t_max * dx, py + t_max * dy);
        return true;
    }

    double getRayLength(const Point& start, const Point& end, int i, int j) const {
        double x1 = static_cast<double>(j);
        double x2 = x1 + 1.0;
        double y1 = static_cast<double>(i);
        double y2 = y1 + 1.0;

        // Ray direction vector
        Point dir = /* TODO */;

        if (std::abs(dir.x) < EPS && std::abs(dir.y) < EPS) {
            return 0.0;
        }

        double t0 = 0.0;
        double t1 = 1.0;

        if (std::abs(dir.x) < EPS) {
            if (start.x < x1 - EPS || start.x > x2 + EPS) {
                return 0.0;
            }
        } else {
            double tx1 = intersectionParamX(x1, start, dir);
            double tx2 = intersectionParamX(x2, start, dir);
            if (tx1 > tx2) {
                std::swap(tx1, tx2);
            }
            t0 = std::max(t0, tx1);
            t1 = std::min(t1, tx2);
        }

        if (t1 <= t0) {
            return 0.0;
        }

        if (std::abs(dir.y) < EPS) {
            if (start.y < y1 - EPS || start.y > y2 + EPS) {
                return 0.0;
            }
        } else {
            double ty1 = intersectionParamY(y1, start, dir);
            double ty2 = intersectionParamY(y2, start, dir);
            if (ty1 > ty2) {
                std::swap(ty1, ty2);
            }
            t0 = std::max(t0, ty1);
            t1 = std::min(t1, ty2);
        }

        if (t1 <= t0) {
            return 0.0;
        }

        double t_enter = std::max(0.0, t0);
        double t_exit  = std::min(1.0, t1);

        if (t_exit <= t_enter) {
            return 0.0;
        }

        double length = /* TODO: length of the ray should be delta(t) * ||dir|| */;

        if (length < 0.0) {
            return 0.0;
        }
        return length;
    }

    double computeProjection(double angle, double offset) const {
        Point start, end;
        if (!buildRaySegment(angle, offset, start, end)) {
            return 0.0;
        }

        double sum = 0.0;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                double len = getRayLength(start, end, i, j);
                sum += /* TODO: see the formula for projection */;
            }
        }
        return sum;
    }

public:
    CTReconstructor(int size) : n(size) {
        matrix.assign(n, std::vector<double>(n, 0.0));
    }

    void readInput() {
        int numAngles;
        std::cin >> numAngles;

        angles.resize(numAngles);
        for (int i = 0; i < numAngles; ++i) {
            double angleDeg;
            std::cin >> angleDeg;
            angles[i] = deg2rad(angleDeg);
        }

        projections.resize(numAngles);
        for (int a = 0; a < numAngles; ++a) {
            int rayCount;
            std::cin >> rayCount;
            projections[a].resize(rayCount);
            for (int r = 0; r < rayCount; ++r) {
                std::cin >> projections[a][r];
            }
        }
    }

    void reconstruct(int iterations) {
        for (int iter = 0; iter < iterations; ++iter) {
            for (std::size_t a = 0; a < angles.size(); ++a) {
                const std::size_t numRays = projections[a].size();

                for (std::size_t r = 0; r < numRays; ++r) {
                    double offset = 0.0;
                    if (numRays > 1) {
                        offset = -static_cast<double>(n) +
                                 2.0 * static_cast<double>(n) * static_cast<double>(r) /
                                 static_cast<double>(numRays - 1);
                    }

                    double measured   = /* TODO */;
                    double calculated = /* TODO */;
                    double error      = measured - calculated;

                    if (std::abs(error) < EPS) {
                        continue;
                    }
                    /* TODO: update the guessed matrix */
                }
            }
        }
    }

private:
    void updateMatrix(double angle, double offset, double error) {
        const double alpha = 0.1;

        Point start, end;
        if (!buildRaySegment(angle, offset, start, end)) {
            return;
        }

        double total_length_sq = 0.0;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                double len = getRayLength(start, end, i, j);
                total_length_sq += len * len;
            }
        }

        if (total_length_sq < EPS) {
            return;
        }

        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                double len = getRayLength(start, end, i, j);
                if (len > EPS) {
                    matrix[i][j] += alpha * error * len / total_length_sq;
                }
            }
        }
    }

public:
    void printMatrix() const {
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                std::cout << std::setw(6) << std::setprecision(6)
                          << std::fixed << matrix[i][j] << " ";
            }
            std::cout << "\n";
        }
    }
};

int main() {
    int n;
    std::cin >> n;

    CTReconstructor ct(n);
    ct.readInput();
    ct.reconstruct(5000);
    ct.printMatrix();

    return 0;
}
中文 English

C++ Online Judge | 开源项目 by Prof. Kaifeng Lyu