VDB Curvature Flow

[Librarian]

/*
 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");
float pc_radius = ch("pc_radius");
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
        vector gradient = volumegradient(2, 0, curr_pos);
        float sample = volumesample(2, 0, curr_pos);
        curr_pos += normalize(-gradient) * sample;
        curr_normal = normalize(-gradient);
        old_pt = new_pt;
        new_pt = addpoint(geoself(), curr_pos);
        
        // 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");
        prim = addprim(geoself(), "polyline");
        addvertex(geoself(), prim, old_pt);
        addvertex(geoself(), prim, new_pt);
}

 

[papsphilip]

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?

 

[konstantin magnus]

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);
vector perp = normalize(cross(grad, {0,1,0}));
vector dir_y = normalize(cross(perp, grad));

 

volume_dir.hiplc

[papsphilip]

1 minute ago, konstantin magnus said:

Hi @papsphilip,

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

vector grad = volumegradient(1, 0, v@P);
vector perp = normalize(cross(grad, {0,1,0}));
vector dir_y = normalize(cross(perp, grad));

 

volume_dir.hiplc

@konstantin magnus thank you!!