# VDB Curvature Flow

## Recommended Posts

```/*
1. Create normal for input points.
2. Cross product of N and a tangent "z-vector",
rotate this by X degrees. Z-vector points in inclination dir.
3. Create a new point by pushing it along the generated direction
for dist_factor.
4. Move the new point to the surface and connect to the previous point.
5.     Create attributes.
6. Repeat.
*/

float dist_factor = ch("dist_factor");
float rotation_degrees = ch("degrees");
float per_step_rotation = ch("per_step_rot");
int steps = chi("steps");
int diagonal = chi("diagonal");

vector find_z_vec(vector pos; float pcrad; int diagonal) {
int handle = pcopen(1, "P", pos, pcrad, 1);
if (pciterate(handle)) {
int found_pt;
vector found_pos;
pcimport(handle, "point.number", found_pt);
pcimport(handle, "P", found_pos);

int neighbors[] = neighbours(1, found_pt);
int n_neighbor_arr[];
int cross_neighbors[];
vector max_y_pos = set(0, -1000, 0);

foreach (int n_pt; neighbors){
vector n_pos = point(1, "P", n_pt);

if (n_pos.y > max_y_pos.y) {
max_y_pos = n_pos;
}
if (diagonal){
int a_neighbors[] = neighbours(1, n_pt);
foreach(int a_pt; a_neighbors){
append(n_neighbor_arr, a_pt);
}
}
}

// Activating diagonal also checks the diagonal neighbors for a maximum y value
if (diagonal) {
foreach(int m_index; int m_pt; n_neighbor_arr){
foreach(int k_index; int k_pt; n_neighbor_arr){
if (k_index > m_index && m_pt == k_pt && m_pt != found_pt){
append(cross_neighbors, m_pt);
}
}
}
foreach(int cross_pt; cross_neighbors){
vector n_pos = point(1, "P", cross_pt);
if (n_pos.y > max_y_pos.y){
max_y_pos = n_pos;
}
}
}

pcclose(handle);
return normalize(found_pos - max_y_pos);

} else {
return set(0, 0, 0);
}
}

int old_pt, new_pt = @ptnum, prim;
vector curr_normal = @N;
vector curr_pos = @P;
for (int i = 0; i < steps; ++i){

// find the tangent
vector zvec = find_z_vec(curr_pos, pc_radius, diagonal);
if (length(zvec) == 0.0) //no max value found => no point on the surface found
break;

// basic translation of the point
vector pre_rot_vec = cross(curr_normal, zvec);
vector4 quat = quaternion(radians(rotation_degrees + per_step_rotation * i), curr_normal);
vector rotated = normalize(qrotate(quat, pre_rot_vec));

curr_pos = curr_pos + dist_factor * rotated;

// this pushes the point to the surface
float sample = volumesample(2, 0, curr_pos);
old_pt = new_pt;

// set attributes and create a spline, grad attribute is used for the color ramp
setpointattrib(geoself(), "N", new_pt, curr_normal, "set");
setpointattrib(geoself(), "grad", new_pt, i / (float) steps, "set");
}
```

##### Share on other sites

is this for calculating the gradient y component of a volume? you are using points so i guess not. its working with a polygon mesh input 0 and a volume in 1. Could you provide an example scene?

##### Share on other sites

Hi @papsphilip,

On 23.9.2021 at 9:58 PM, papsphilip said:

gradient y component of a volume

this would be the double cross product of the volume's gradient:

```vector grad = volumegradient(1, 0, v@P);

volume_dir.hiplc

• 1

##### Share on other sites
1 minute ago, konstantin magnus said:

Hi @papsphilip,

this would be the double cross product of the volume's gradient:

```

@konstantin magnus thank you!!

## Create an account

Register a new account