Jump to content

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");
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);
}

 

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

 

volume_dir.hiplc

  • Thanks 1
Link to comment
Share on other sites

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!!

Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

×
×
  • Create New...