Function generate_image

examples/pde/bhrt.cpp:917–1066

* @brief Returns an rgb image containing the raytraced black hole from the * camera rays, spacetime metric, objects living in the space, and background * * @param initial_pos initial position from where the rays are launched * @param initial_vel initial velocities (directions) the rays have * @param objects the objects the rays can collide with * @param background the background of the scene

 * @return af::array
 */
af::array generate_image(const af::array& initial_pos,
                         const af::array& initial_vel,
                         const std::vector<std::unique_ptr<Object> >& objects,
                         const Background& background, uint32_t width,
                         uint32_t height, double time, double tol,
                         uint32_t checks = 10) {
    uint32_t lines = initial_pos.dims()[1];

    auto def_step = 0.5 * pow(tol, 0.25);
    auto dt       = af::constant(def_step, 1, lines, f64);
    auto t        = af::constant(0.0, 1, lines, f64);
    auto index    = af::iota(lines);
    auto selected = t < time;

    auto result = af::constant(0, lines, 3, f32);

    auto pos = initial_pos;
    auto vel = initial_vel;

    af::Window window{(int)width, (int)height, "Black Hole Raytracing"};

    af::array bg_col = af::constant(0.f, lines, 3);
    af::array begin_pos, end_pos;
    af::array bh_nohit;

    if (scene != Scene::ROTATE_BH)
        begin_pos = sph_to_cart_position(pos(af::seq(1, 3), af::span));
    else
        begin_pos = oblate_to_cart_position(pos(af::seq(1, 3), af::span));
    end_pos = begin_pos;

    int i = 0;

    while (t.dims()[1] != 0 && af::anyTrue<bool>(t < time) &&
           af::anyTrue<bool>(dt != 0.0)) {
        // Displays the current progress and approximate time needed to finish
        // it
        status_bar((lines - t.dims()[1]) * time +
                       af::sum<double>(af::clamp(t, 0.0, time)),
                   time * lines, "Progress:");

        // RK34 method for second order differential equation
        auto dt2 = dt * dt;
        auto k1  = geodesics(pos, vel);
        auto k2  = geodesics(pos + vel * dt / 4.0 + k1 * dt2 / 32.0,
                             vel + k1 * dt / 4.0);
        auto k3  = geodesics(pos + vel * dt / 2.0 + (k1 + k2) * dt2 / 16.0,
                             vel + k2 * dt / 2.0);
        auto k4  = geodesics(pos + vel * dt + (k1 - k2 + k3 * 2.0) * dt2 / 4.0,
                             vel + (k1 - k2 * 2.0 + 2.0 * k3) * dt);

        auto diff4 = (k1 + k2 * 8.0 + k3 * 2.0 + k4) / 24.0;
        auto diff3 = (k2 * 8.0 + k4) / 18.0;

        auto err    = (af::max)(af::abs(diff4 - diff3), 0) * dt2;
        auto maxerr = tol * (1.0 + (af::max)(af::abs(pos), 0));

        auto rdt = af::constant(0, 1, dt.dims()[1], f64);