Generate vtk file of fractal.

genVtk.cpp

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <vector>
#include <cmath>
#include "nanort.h"

typedef nanort::real3<int> int3;
typedef nanort::real3<float> float3;
typedef nanort::real3<double> double3;

typedef struct {
    float3 center;
    float r;
} Sphere;

typedef struct {
    float3 normal;
    float3 origin;
} Plane;

template<typename T>
inline T vdistance(nanort::real3<T> a, nanort::real3<T> b) {
    nanort::real3<T> d = a - b;
    return std::sqrt(nanort::vdot(d, d));
}

inline void sphereInvert(float3 &pos, float &dr, const Sphere s) {
    float3 diff = pos - s.center;
    float lenSq = nanort::vdot(diff, diff);
    float k = (s.r * s.r) / lenSq;
    dr *= k; // (r * r) / lenSq
    pos = diff * k + s.center;
}

inline float distSphere(const float3 pos, const Sphere s) {
    return vdistance(pos, s.center) - s.r;
}

inline float distPlane(const float3 pos, const Plane p) {
    return vdot(pos - p.origin, p.normal);
}

inline float distInfSphairahedron(const float3 pos,
                                  const std::vector<Sphere> spheres,
                                  const std::vector<Plane> planes,
                                  const Plane divPlane) {
    float d = -1;
    d = std::max(distPlane(pos, planes[0]), d);
    d = std::max(distPlane(pos, planes[1]), d);
    d = std::max(distPlane(pos, planes[2]), d);
    d = std::max(distPlane(pos, divPlane), d);

    d = std::max(-distSphere(pos, spheres[0]), d);
    d = std::max(-distSphere(pos, spheres[1]), d);
    d = std::max(-distSphere(pos, spheres[2]), d);
    return d;
}

const float3 tP1(1 + 0.5, 0, -std::sqrt(3.) + std::sqrt(3.) * 0.5);
const float3 tP2(-(1 + 0.5), 0, -std::sqrt(3.) + std::sqrt(3.) * 0.5);
const float3 tP3(-0.5 + 0.5, 0, std::sqrt(3.) * 0.5 + std::sqrt(3.) * 0.5);
const Plane testP1 = {float3(1, 0, 0),
                      tP1 };
//vnormalize((tP1 + float3(0, 0, -std::sqrt(3.))) * 0.5)
const Plane testP2 = {float3(0.5,
                             0,
                             -0.8660254037844388),
                      tP1 };
const Plane testP3 = {float3(-0.5,
                             0,
                             -0.8660254037844388),
                      tP2 };
const Plane testP4 = {float3(-1, 0, 0),
                      tP2 };
const Plane testP5 = {float3(0.5,
                             0,
                             0.8660254037844388),
                      tP3 };
const Plane testP6 = {float3(-0.5,
                             0,
                             0.8660254037844388),
                      tP3 };

std::vector<Plane> testPlanes = { testP1, testP2, testP3, testP4, testP5, testP6 };

const int MAX_ITER_COUNT = 1000;
float distFunc(float3 pos,
               const std::vector<Sphere> spheres,
               const std::vector<Plane> planes,
               const Plane divPlane) {
    bool outside = false;
    std::for_each(testPlanes.begin(), testPlanes.end(),
                  [&](Plane p){
                      pos = pos - p.origin;
                      float d = vdot(pos, p.normal);
                      if (d > 0.) {
                          outside = true;
                      }
                      pos = pos + p.origin;
                  });
    if(outside) return 10.;
    int invNum = 0;
    float dr = 1.0;
    for(int n = 0; n < MAX_ITER_COUNT; n++) {
        bool inFund = true;
        std::for_each(spheres.begin(), spheres.end(),
                      [&](Sphere s){
                          if (vdistance(pos, s.center) < s.r) {
                              //std::cout << pos.x() << std::endl;
                              sphereInvert(pos, dr, s);
                              //std::cout << pos.x() << std::endl;
                              invNum++;
                              inFund = false;
                          }
                      });
        std::for_each(planes.begin(), planes.end(),
                      [&](Plane p){
                          pos = pos - p.origin;
                          float d = vdot(pos, p.normal);
                          if (d > 0.) {
                              pos = pos - 2.f * d * p.normal;
                              invNum++;
                              inFund = false;
                          }
                          pos = pos + p.origin;
                      });
        if (inFund) break;
    }
    return distInfSphairahedron(pos, spheres, planes, divPlane) / dr;
}

int main() {
    std::ofstream ofs("testIIS.vtk");
    if(!ofs) {
        std::cerr << "error" << std::endl;
		std::exit(EXIT_FAILURE);
    }

    float3 offset = float3(0.5, 0, std::sqrt(3.) * 0.5);
//    offset = float3(0, 0, 0);
    Sphere s1 = { float3(0.25450630091522675,
                         0,
                         0),
                  0.7454936990847733 };
    s1.center = s1.center + offset;
    Sphere s2 = { float3(-0.3116633792528053,
                         0.6053931133878944,
                         0.5398168077244666),
                  0.37667324149438935 };
    s2.center = s2.center + offset;
    Sphere s3 = { float3(-0.12608782704164367,
                         1.2165336555165491,
                         -0.21839052265208383),
                  0.7478243459167127 };
    s3.center = s3.center + offset;
    std::vector<Sphere> prismSpheres = { s1, s2, s3 };

    Plane p1 = { float3(0.5,
                        0,
                        0.8660254037844388),
                 float3(0,
                        5,
                        0.5773502691896258) };
    p1.origin = p1.origin + offset;
    Plane p2 = { float3(0.5,
                        0,
                        -0.8660254037844388),
                 float3(0,
                        3,
                        -0.5773502691896258) };
    p2.origin = p2.origin + offset;
    Plane p3 = { float3(-1, 0, 0),
                 float3(-0.5, 0, 1) };
    p3.origin = p3.origin + offset;
    std::vector<Plane> prismPlanes = { p1, p2, p3 };

    Plane divPlane = { float3(0.4969732028017572,
                              0.8183202716219945,
                              0.2887378893554992),
                       float3(0.9999999999999948,
                              -1.3100631690576847e-14,
                              7.91033905045424e-15) };
    divPlane.origin = divPlane.origin + offset;

    float3 origin = float3(-1.8, -1.0, -1.9);
    float3 n = float3(1.8, 1.3, 1.9);
    float3 s = float3(0.05, 0.05, 0.05);
    int3 dim = int3(int((n.x() - origin.x()) / s.x()),
                    int((n.y() - origin.y()) / s.y()),
                    int((n.z() - origin.z()) / s.z()));
    int numPoints = dim.x() * dim.y() * dim.z();

    std::vector<int> invNums;
    invNums.resize(numPoints);
    std::vector<float> distances;
    distances.resize(numPoints);
    int i = 0;
    std::cout << "dim "
              << dim.x() << " "
              << dim.y() << " "
              << dim.z() << std::endl;
    std::cout << "numPoints " << numPoints << std::endl;

    std::cout << "start" << std::endl;
    for(int zi = 0; zi < dim.z(); zi++) {
        std::cout << zi << std::endl;
        for(int yi = 0; yi < dim.y(); yi++) {
            for(int xi = 0; xi < dim.x(); xi++) {
                if (yi == 0) {
                    distances[i] = 0.;
                    i++;
                    continue;
                }
                float3 p = float3(origin.x() + s.x() * float(xi),
                                  origin.y() + s.y() * float(yi),
                                  origin.z() + s.z() * float(zi));
                float d = distFunc(p, prismSpheres, prismPlanes, divPlane);
                if (d > 0.001) {
                    distances[i] = 0.;
                } else {
                    distances[i] = 10.;
                }
                i++;
            }
        }
    }
    std::cout << "done" << std::endl;

    ofs << "# vtk DataFile Version 3.0" << std::endl;
    ofs << "vtk output" << std::endl;
    ofs << "ASCII" << std::endl;
    ofs << "DATASET STRUCTURED_POINTS" << std::endl;
    ofs << "DIMENSIONS "
        << int((n.x() - origin.x()) / s.x()) << " "
        << int((n.y() - origin.y()) / s.y()) << " "
        << int((n.z() - origin.z()) / s.z())
        << std::endl;
    ofs << "ORIGIN "
        << origin.x() << " "
        << origin.y() << " "
        << origin.z() << std::endl;
    ofs << "SPACING "
        << s.x() << " "
        << s.y() << " "
        << s.z() << std::endl;
    ofs << "POINT_DATA "
        << numPoints
        << std::endl;
    ofs << "SCALARS scalars float" << std::endl;
    ofs << "LOOKUP_TABLE default" << std::endl;
    std::cout << "Writing points..." << std::endl;
    std::for_each(distances.begin(), distances.end(),
                  [&](float p) {
                      ofs << p << std::endl;
                  });
    std::cout << "Done" << std::endl;

    ofs.close();
    std::cout << "end" << std::endl;
}