Skip to content

Instantly share code, notes, and snippets.

@soma-arc soma-arc/genVtk.cpp
Created Dec 23, 2017

Embed
What would you like to do?
Generate vtk file of fractal.
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.