examples/cpu/volume_rendering/volume.cpp - main - redgpu/ispc

examples/cpu/volume_rendering/volume.cpp
195 строк6 KB
Dmitry Babokin
01 апр 2023, 03:48
01 апр 2023, 03:4833bd2b4
/*
  Copyright (c) 2011-2023, Intel Corporation

  SPDX-License-Identifier: BSD-3-Clause
*/

#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#define NOMINMAX
#pragma warning(disable : 4244)
#pragma warning(disable : 4305)
#endif

#include "../../common/timing.h"
#include "volume_ispc.h"
#include <algorithm>
#include <cstdlib>
#include <stdio.h>
using namespace ispc;

extern void volume_serial(float density[], int nVoxels[3], const float raster2camera[4][4],
                          const float camera2world[4][4], int width, int height, float image[]);

/* Write a PPM image file with the image */
static void writePPM(float *buf, int width, int height, const char *fn) {
    FILE *fp = fopen(fn, "wb");
    if (!fp) {
        printf("Couldn't open a file '%s'\n", fn);
        exit(1);
    }
    fprintf(fp, "P6\n");
    fprintf(fp, "%d %d\n", width, height);
    fprintf(fp, "255\n");
    for (int i = 0; i < width * height; ++i) {
        float v = buf[i] * 255.f;
        if (v < 0.f)
            v = 0.f;
        else if (v > 255.f)
            v = 255.f;
        unsigned char c = (unsigned char)v;
        for (int j = 0; j < 3; ++j)
            fputc(c, fp);
    }
    fclose(fp);
    printf("Wrote image file %s\n", fn);
}

/* Load image and viewing parameters from a camera data file.
   FIXME: we should add support to be able to specify viewing parameters
   in the program here directly. */
static void loadCamera(const char *fn, int *width, int *height, float raster2camera[4][4], float camera2world[4][4]) {
    FILE *f = fopen(fn, "r");
    if (!f) {
        perror(fn);
        exit(1);
    }
    if (fscanf(f, "%d %d", width, height) != 2) {
        fprintf(stderr, "Unexpected end of file in camera file\n");
        exit(1);
    }

    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            if (fscanf(f, "%f", &raster2camera[i][j]) != 1) {
                fprintf(stderr, "Unexpected end of file in camera file\n");
                exit(1);
            }
        }
    }
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            if (fscanf(f, "%f", &camera2world[i][j]) != 1) {
                fprintf(stderr, "Unexpected end of file in camera file\n");
                exit(1);
            }
        }
    }
    fclose(f);
}

/* Load a volume density file.  Expects the number of x, y, and z samples
   as the first three values (as integer strings), then x*y*z
   floating-point values (also as strings) to give the densities.  */
static float *loadVolume(const char *fn, int n[3]) {
    FILE *f = fopen(fn, "r");
    if (!f) {
        perror(fn);
        exit(1);
    }

    if (fscanf(f, "%d %d %d", &n[0], &n[1], &n[2]) != 3) {
        fprintf(stderr, "Couldn't find resolution at start of density file\n");
        fclose(f);
        exit(1);
    }

    int count = n[0] * n[1] * n[2];
    float *v = new float[count];
    for (int i = 0; i < count; ++i) {
        if (fscanf(f, "%f", &v[i]) != 1) {
            fprintf(stderr, "Unexpected end of file at %d'th density value\n", i);
            fclose(f);
            exit(1);
        }
    }

    fclose(f);
    return v;
}

int main(int argc, char *argv[]) {
    static unsigned int test_iterations[] = {3, 7, 1};
    if (argc < 3) {
        fprintf(stderr, "usage: volume <camera.dat> <volume_density.vol> [ispc iterations] [tasks iterations] [serial "
                        "iterations]\n");
        return 1;
    }
    if (argc == 6) {
        for (int i = 0; i < 3; i++) {
            test_iterations[i] = atoi(argv[3 + i]);
        }
    }

    //
    // Load viewing data and the volume density data
    //
    int width, height;
    float raster2camera[4][4], camera2world[4][4];
    loadCamera(argv[1], &width, &height, raster2camera, camera2world);
    float *image = new float[width * height];

    int n[3];
    float *density = loadVolume(argv[2], n);

    //
    // Compute the image using the ispc implementation; report the minimum
    // time of three runs.
    //
    double minISPC = 1e30;
    for (unsigned int i = 0; i < test_iterations[0]; ++i) {
        reset_and_start_timer();
        volume_ispc(density, n, raster2camera, camera2world, width, height, image);
        double dt = get_elapsed_mcycles();
        printf("@time of ISPC run:\t\t\t[%.3f] million cycles\n", dt);
        minISPC = std::min(minISPC, dt);
    }

    printf("[volume ispc 1 core]:\t\t[%.3f] million cycles\n", minISPC);
    writePPM(image, width, height, "volume-ispc-1core.ppm");

    // Clear out the buffer
    for (int i = 0; i < width * height; ++i)
        image[i] = 0.;

    //
    // Compute the image using the ispc implementation that also uses
    // tasks; report the minimum time of three runs.
    //
    double minISPCtasks = 1e30;
    for (unsigned int i = 0; i < test_iterations[1]; ++i) {
        reset_and_start_timer();
        volume_ispc_tasks(density, n, raster2camera, camera2world, width, height, image);
        double dt = get_elapsed_mcycles();
        printf("@time of ISPC + TASKS run:\t\t\t[%.3f] million cycles\n", dt);
        minISPCtasks = std::min(minISPCtasks, dt);
    }

    printf("[volume ispc + tasks]:\t\t[%.3f] million cycles\n", minISPCtasks);
    writePPM(image, width, height, "volume-ispc-tasks.ppm");

    // Clear out the buffer
    for (int i = 0; i < width * height; ++i)
        image[i] = 0.;

    //
    // And run the serial implementation 3 times, again reporting the
    // minimum time.
    //
    double minSerial = 1e30;
    for (unsigned int i = 0; i < test_iterations[2]; ++i) {
        reset_and_start_timer();
        volume_serial(density, n, raster2camera, camera2world, width, height, image);
        double dt = get_elapsed_mcycles();
        printf("@time of serial run:\t\t\t[%.3f] million cycles\n", dt);
        minSerial = std::min(minSerial, dt);
    }

    printf("[volume serial]:\t\t[%.3f] million cycles\n", minSerial);
    writePPM(image, width, height, "volume-serial.ppm");

    printf("\t\t\t\t(%.2fx speedup from ISPC, %.2fx speedup from ISPC + tasks)\n", minSerial / minISPC,
           minSerial / minISPCtasks);

    return 0;
}