sandbox/Antoonvh/tree/treegen.h
Generate a fractal tree
Stem parameters:
coord stem_pos;
double stem_rad = 1.; //Base size
Length-over-radius ratio, and number of fractal levels. The stem is at level = 0.
double LoR = 20., LoR_spread = 5.; //LoR_i = LoR + noise()*LoR_spread
int levels = 1;
Branch parameters: Number of branches will be rounded to the nearest integer after the spreading randomization. The new_angle
parameter is the angle relative to the parent branch (100% randomized rotation).
double fractal_factor = 0.7, fractal_factor_spread = 0.1;
double split_branches = 2.1, split_branches_spread = 0.5; //Starting at the end of parent
double side_branches = 4.1, side_branches_spread = 0.5; //Starting along the parent
double side_start = 0.65, side_start_spread = .3; //Relative along-parent position
double new_angle = 40.*pi/180., new_angle_spread = 10.*pi/180.;
That concludes the inputs.
typedef struct {
coord start, end;
double R;
int parent;
} Branch;
unsigned int nb = 0; //Number of segments
coord cross (coord a, coord b) {
coord new;
new.x = a.y*b.z - a.z*b.y;
new.y = a.z*b.x - a.x*b.z;
new.z = a.x*b.y - a.y*b.x;
return new;
}
void normalize3 (coord * n)
{
double norm = 0.;
foreach_dimension(3)
norm += sq(n->x);
norm = sqrt(norm);
foreach_dimension(3)
n->x /= norm;
}
unsigned int add_bones (Branch * branches, unsigned int p, unsigned int j) {
Branch par = branches[p];
int indp = j;
int sidenr = (int)round(side_branches + side_branches_spread*noise()) + 0.5;
double rndm_ngl = pi*noise();
coord randomvec = {noise(), noise(), noise()};
for (int m = 0; m < sidenr; m++) {
j++;
branches[j].parent = indp;
double sd_strt = side_start + noise()*side_start_spread;
foreach_dimension(3)
branches[j].start.x = par.start.x + sd_strt*(par.end.x-par.start.x);
double frcl_fctr = fractal_factor + fractal_factor_spread*noise();
branches[j].R = par.R*frcl_fctr;
double length = (LoR + noise()*LoR_spread)*branches[j].R;
double ngl = new_angle + new_angle_spread*noise();
rndm_ngl += 2.*pi/((double)sidenr) + new_angle_spread*noise();
coord vec;
foreach_dimension(3) {
vec.x = (par.end.x - par.start.x);
}
We need to rotate this vector with angle ngl
and rndm_ngl
(pitch and yaw)
This was helpfull:
tparker, Rotate a vector by a randomly oriented angle, URL (version: 2017-09-30): https://scicomp.stackexchange.com/q/27969
normalize3 (&vec);
coord vecx = cross (vec, randomvec);
normalize3 (&vecx);
coord vecy = cross (vec, vecx);
coord vecn;
foreach_dimension(3)
vecn.x = (sin(ngl)*cos(rndm_ngl)*vecx.x +
sin(ngl)*sin(rndm_ngl)*vecy.x +
cos(ngl)*vec.x);
foreach_dimension(3)
branches[j].end.x = branches[j].start.x + vecn.x*length;
foreach_dimension(3)
branches[j].start.x += vecn.x*par.R; //away from parent center line
}
int splitnr = (int)round(split_branches + split_branches_spread*noise()) + 0.5;
rndm_ngl = pi*noise();
for (int m = 0; m < splitnr; m++) {
j++;
branches[j].parent = indp;
branches[j].start = par.end;
double frcl_fctr = fractal_factor + fractal_factor_spread*noise();
branches[j].R = par.R*frcl_fctr;
double length = (LoR + noise()*LoR_spread)*branches[j].R;
double ngl = new_angle + new_angle_spread*noise();
rndm_ngl += 2*pi/((double)splitnr) + new_angle_spread*noise();
coord vec;
foreach_dimension(3)
vec.x = (par.end.x - par.start.x);
normalize3 (&vec);
coord vecx = cross (vec, randomvec);
normalize3 (&vecx);
coord vecy = cross (vec, vecx);
coord vecn;
foreach_dimension(3)
vecn.x = (sin(ngl)*cos(rndm_ngl)*vecx.x +
cos(ngl)*sin(rndm_ngl)*vecy.x +
cos(ngl)*vec.x);
foreach_dimension(3)
branches[j].end.x = branches[j].start.x + vecn.x*length;
foreach_dimension(3)
branches[j].start.x += vecn.x*par.R; //away from parent center line
}
return splitnr + sidenr;
}
Branch * tree_skeleton () {
Branch * branches;
//array or branches
int maxb = (round(split_branches + split_branches_spread) +
round(side_branches + side_branches_spread));
unsigned int maxnl = (unsigned int)(pow(maxb, levels + 1) + 3);
branches = (Branch*) malloc(sizeof(Branch)*maxnl);
// init stem
int j = 0;
branches[j].start = stem_pos;
branches[j].R = stem_rad;
coord ey = {0, 1, 0};
double length = (LoR + noise()*LoR_spread)*branches[j].R;
foreach_dimension(3)
branches[j].end.x = stem_pos.x + length*ey.x;
branches[j].parent = 0;
nb ++;
int pc = 0;
for (int l = 0 ; l < levels; l++) {
// At this level, all branches are parents
while (pc < nb) {
int n = add_bones (branches, pc, j); //parent pc to child j
pc++;
j += n;
}
nb = j + 1;
}
return branches;
}
Skeleton to volume fractions
#include "PointTriangle.h" //PointSegmentDistance exists!
#include "fractions.h"
struct tree_int {
Branch * branches;
scalar c;
face vector fs;
scalar J;
double smooth;
scalar * alist; //For adaptation
double * crit; //...
int minlevel;
int maxlevel;
scalar * ulist;
int stop;
double stopc;
};
#define SMOOTHR ((d - branches[j].R) > 0 ?exp (-sq((d - branches[j].R)/(branches[j].R*ti.smooth))) : 1)
//#define SMOOTHR ((d - branches[j].R) > 0 ? ((d - branches[j].R)/branches[j].R < pi/2 ? (cos (d - branches[j].R)/branches[j].R) : 0) : 1)
trace
void tree_interface (struct tree_int ti) {
Branch * branches = ti.branches;
vertex scalar phi[], phiR[];
scalar J = automatic (ti.J);
if (!ti.smooth) {
foreach_vertex() {
phi[] = HUGE;
for (int j = 0; j < nb; j++) {
coord pos = (coord){x, y, z};
coord p0 = branches[j].start;
coord p1 = branches[j].end;
coord Sc;
double sP;
double d = sqrt(PointSegmentDistance (&pos, &p0, &p1, &Sc, &sP))
- branches[j].R;
if (d < phi[]) {
phi[] = d;
if (ti.J.i)
J[] = j;
}
}
}
} else {
foreach_vertex() {
phi[] = 0;
double mind = HUGE, phiR = 0;
for (int j = 0; j < nb; j++) {
coord pos = (coord){x, y, z};
coord p0 = branches[j].start;
coord p1 = branches[j].end;
coord Sc;
double sP;
double d = sqrt(PointSegmentDistance (&pos, &p0, &p1, &Sc, &sP));
phi[] += d > 0 ? sq(1/d) : HUGE ;
phiR += SMOOTHR*branches[j].R;
if (ti.J.i)
if (d < mind) {
mind = d;
J[] = j;
}
}
phi[] = sqrt(1./phi[]) - phiR;
}
}
boundary ({phi});
fractions (phi, ti.c, ti.fs);
}
Skeleton to fractions with refinement
void refine_vertex (Point point, scalar s) {
for (int i = 0; i < 2; i++)
for (int j = 0; j < 1 + (dimension > 1); j++)
for (int k = 0; k < 1 + (dimension > 2); k++) {
fine (s, 2*i, 2*j, 2*k) = s[i, j, k];
fine (s, 2*i + 1, 2*j , 2*k) = nodata; //1
#if (dimension > 1)
fine (s, 2*i , 2*j + 1, 2*k) = nodata; //1
fine (s, 2*i + 1, 2*j + 1, 2*k) = nodata; //2
#if (dimension > 2)
fine (s, 2*i , 2*j , 2*k + 1) = nodata; //1
fine (s, 2*i + 1, 2*j , 2*k + 1) = nodata; //2
fine (s, 2*i , 2*j + 1, 2*k + 1) = nodata; //2
fine (s, 2*i +1 , 2*j + 1, 2*k + 1) = nodata; //3
#endif
#endif
}
}
void tree_interface_adapt (struct tree_int ti) {
Branch * branches = ti.branches;
vertex scalar phi[];
phi.refine = refine_vertex;
foreach_vertex()
phi[] = nodata;
scalar * update = NULL;
if (!ti.ulist)
update = list_copy (all);
else {
update = list_copy (ti.ulist);
update = list_add(update, phi);
}
int it = 0, computed;
do {
computed = 0;
if (ti.stopc)
ti.stop = grid->tn/(int)ti.stopc;
if (!ti.smooth) {
foreach_vertex(reduction (+:computed)) {
if (phi[] == nodata) {
computed++;
for (int j = 0; j < nb; j++) {
coord pos = (coord){x, y, z};
coord p0 = branches[j].start;
coord p1 = branches[j].end;
coord Sc;
double sP;
double d = sqrt(PointSegmentDistance (&pos, &p0, &p1, &Sc, &sP))
- branches[j].R;
if (d < phi[])
phi[] = d;
}
}
}
} else {
foreach_vertex(reduction (+:computed)) {
if (phi[] == nodata) {
computed++;
phi[] = 0;
double phiR = 0;
for (int j = 0; j < nb; j++) {
coord pos = (coord){x, y, z};
coord p0 = branches[j].start;
coord p1 = branches[j].end;
coord Sc;
double sP;
double d = sqrt(PointSegmentDistance (&pos, &p0, &p1, &Sc, &sP));
phi[] += d > 0 ? sq(1/d) : HUGE;
phiR += SMOOTHR*branches[j].R;
}
phi[] = sqrt(1./phi[]) - phiR;
}
}
}
boundary ({phi});
fractions (phi, ti.c, ti.fs);
scalar c = ti.c;
face vector fss = ti.fs;
boundary ({c, fss});
#if EMBED
if (fractions_cleanup (ti.c, ti.fs))
boundary ({c, fss});
#endif
if (pid() == 0)
printf ("# Tree Reconstruction Adapt Iteration %d: %ld cells\n"
"# and %d newly computed vertices\n", ++it, grid->tn, computed);
} while (adapt_wavelet (ti.alist, ti.crit, ti.maxlevel, ti.minlevel, update).nf > ti.stop);
scalar c = ti.c;
face vector fss = ti.fs;
boundary ({c, fss});
// Last time for R
if (ti.J.i) {
tree_interface (ti);
boundary ({c, fss});
#if EMBED
if (fractions_cleanup (ti.c, ti.fs))
boundary ({c, fss});
#endif
}
}