@@ -174,6 +174,9 @@ def _simple_grid_3d_octree_regular():
174174 data_shape = data_shape ,
175175 options = options ,
176176 )
177+
178+ intersection_xyz , valid_edges = _grab_xyz_edges (octree_level_for_surface )
179+
177180 dc_meshes : List [DualContouringMesh ] = compute_dual_contouring (dc_data )
178181 dc_data = dc_meshes [0 ].dc_data
179182 valid_voxels = dc_data .valid_voxels
@@ -321,11 +324,11 @@ def test_compute_dual_contouring_several_meshes(simple_model_3_layers, simple_gr
321324
322325
323326@pytest .mark .skipif (BackendTensor .engine_backend != AvailableBackends .numpy , reason = "Only numpy supported" )
324- def test_find_edges_intersection_step_by_step (simple_model , simple_grid_3d_octree ):
327+ def test_find_edges_intersection_pro (simple_model , simple_grid_3d_octree ):
325328 # region Test find_intersection_on_edge
326329 spi , ori_i , options , data_shape = simple_model
327330 ids = np .array ([1 , 2 ])
328- grid_0_centers = copy . deepcopy ( simple_grid_3d_octree )
331+ grid_0_centers = simple_grid_3d_octree
329332 interpolation_input = InterpolationInput (spi , ori_i , grid_0_centers , ids )
330333
331334 options .number_octree_levels = 5
@@ -437,111 +440,6 @@ def test_find_edges_intersection_step_by_step(simple_model, simple_grid_3d_octre
437440 plot_label = False , plot_marching_cubes = False )
438441
439442
440- @pytest .mark .skipif (BackendTensor .engine_backend != AvailableBackends .numpy , reason = "Only numpy supported" )
441- def test_find_edges_intersection_pro (simple_model , simple_grid_3d_octree ):
442- # region Test find_intersection_on_edge
443- spi , ori_i , options , data_shape = simple_model
444- ids = np .array ([1 , 2 ])
445- grid_0_centers = simple_grid_3d_octree
446- interpolation_input = InterpolationInput (spi , ori_i , grid_0_centers , ids )
447-
448- options .compute_scalar_gradient = True
449- octree_list = interpolate_n_octree_levels (interpolation_input , options , data_shape )
450-
451- last_octree_level : OctreeLevel = octree_list [- 1 ]
452-
453- sfsp = last_octree_level .last_output_center .scalar_field_at_sp
454- # sfsp = np.append(sfsp, -0.1)
455- xyz_on_edge , valid_edges = find_intersection_on_edge (last_octree_level .grid_centers .corners_grid .values , last_octree_level .last_output_center .exported_fields .scalar_field [last_octree_level .last_output_center .grid .corners_grid_slice ], sfsp , )
456- # endregion
457-
458- # region Get Normals
459- interpolation_input .set_temp_grid (EngineGrid .from_xyz_coords (xyz_on_edge ))
460-
461- output_on_edges = interpolate_and_segment (interpolation_input , options , data_shape .tensors_structure )
462- # stack gradients output_on_edges.exported_fields.gx_field
463- gradients = np .stack (
464- (output_on_edges .exported_fields .gx_field ,
465- output_on_edges .exported_fields .gy_field ,
466- output_on_edges .exported_fields .gz_field ), axis = 0 ).T
467-
468- # endregion
469-
470- # region Prepare data for vectorized QEF
471-
472- n_edges = valid_edges .shape [0 ]
473-
474- # Coordinates for all posible edges (12) and 3 dummy normals in the center
475- xyz = np .zeros ((n_edges , 15 , 3 ))
476- normals = np .zeros ((n_edges , 15 , 3 ))
477-
478- xyz [:, :12 ][valid_edges ] = xyz_on_edge
479- normals [:, :12 ][valid_edges ] = gradients
480-
481- BIAS_STRENGTH = 0.1
482-
483- xyz_aux = np .copy (xyz [:, :12 ])
484-
485- # Numpy zero values to nans
486- xyz_aux [np .isclose (xyz_aux , 0 )] = np .nan
487- # Mean ignoring nans
488- mass_points = np .nanmean (xyz_aux , axis = 1 )
489-
490- xyz [:, 12 ] = mass_points
491- xyz [:, 13 ] = mass_points
492- xyz [:, 14 ] = mass_points
493-
494- normals [:, 12 ] = np .array ([BIAS_STRENGTH , 0 , 0 ])
495- normals [:, 13 ] = np .array ([0 , BIAS_STRENGTH , 0 ])
496- normals [:, 14 ] = np .array ([0 , 0 , BIAS_STRENGTH ])
497-
498- # Remove unused voxels
499- bo = valid_edges .sum (axis = 1 , dtype = bool )
500- xyz = xyz [bo ]
501- normals = normals [bo ]
502-
503- # Compute LSTSQS in all voxels at the same time
504- A1 = normals
505- b1 = xyz
506- bb1 = (A1 * b1 ).sum (axis = 2 )
507- s1 = np .einsum ("ijk, ilj->ikl" , A1 , np .transpose (A1 , (0 , 2 , 1 )))
508- s2 = np .linalg .inv (s1 )
509- s3 = np .einsum ("ijk,ik->ij" , np .transpose (A1 , (0 , 2 , 1 )), bb1 )
510- v_pro = np .einsum ("ijk, ij->ik" , s2 , s3 )
511-
512- # endregion
513-
514- # endregion
515-
516- # region triangulate
517- grid_centers = last_octree_level .grid_centers
518- valid_voxels = valid_edges .sum (axis = 1 , dtype = bool )
519-
520- temp_ids = octree_list [- 1 ].last_output_center .ids_block # ! I need this because setters in python sucks
521- temp_ids [valid_voxels ] = 5
522- octree_list [- 1 ].last_output_center .ids_block = temp_ids # paint valid voxels
523-
524- dc_data = DualContouringData (
525- xyz_on_edge = xyz_on_edge ,
526- xyz_on_centers = grid_centers .values ,
527- dxdydz = grid_centers .octree_dxdydz ,
528- valid_edges = valid_edges ,
529- exported_fields_on_edges = None ,
530- n_surfaces_to_export = data_shape .tensors_structure .n_surfaces
531- )
532-
533- indices = triangulate_dual_contouring (dc_data )
534- # endregion
535-
536- if plot_pyvista or False :
537- # ! I leave this test for the assert as comparison to the other implementation. The model looks bad
538- # ! with this level of BIAS
539- center_mass = xyz [:, 12 :].reshape (- 1 , 3 )
540- normals = normals [:, 12 :].reshape (- 1 , 3 )
541- _plot_pyvista (last_octree_level , octree_list , simple_model , ids , grid_0_centers ,
542- xyz_on_edge , gradients , a = center_mass , b = normals ,
543- v_pro = v_pro , indices = indices , plot_marching_cubes = True
544- )
545443
546444
547445@pytest .mark .skipif (BackendTensor .engine_backend != AvailableBackends .numpy , reason = "Only numpy supported" )
0 commit comments