@@ -418,6 +418,44 @@ void Cell::refine_func(node_map_t& nodes, function test_func, double *xs, double
418418 }
419419}
420420
421+ void Cell::refine_image (node_map_t & nodes, double * image, int_t *shape_cells, double *xs, double *ys, double *zs, bool diag_balance){
422+ // early exit if my level is higher than or equal to target
423+ if (level == max_level){
424+ return ;
425+ }
426+ int_t start_ix = points[0 ]->location_ind [0 ]/2 ;
427+ int_t start_iy = points[0 ]->location_ind [1 ]/2 ;
428+ int_t start_iz = n_dim == 2 ? 0 : points[0 ]->location_ind [2 ]/2 ;
429+ int_t end_ix = points[3 ]->location_ind [0 ]/2 ;
430+ int_t end_iy = points[3 ]->location_ind [1 ]/2 ;
431+ int_t end_iz = n_dim == 2 ? 1 : points[7 ]->location_ind [2 ]/2 ;
432+ int_t nx = shape_cells[0 ];
433+ int_t ny = shape_cells[1 ];
434+ int_t nz = shape_cells[2 ];
435+
436+ int_t i_image = (nx * ny) * start_iz + nx * start_iy + start_ix;
437+ double val_start = image[i_image];
438+ bool all_unique = true ;
439+
440+ // if any of the image data contained in the cell are different, subdivide myself
441+ for (int_t iz=start_iz; iz<end_iz && all_unique; ++iz)
442+ for (int_t iy=start_iy; iy<end_iy && all_unique; ++iy)
443+ for (int_t ix=start_ix; ix<end_ix && all_unique; ++ix){
444+ i_image = (nx * ny) * iz + nx * iy + ix;
445+ all_unique = image[i_image] == val_start;
446+ }
447+
448+ if (!all_unique){
449+ if (is_leaf ()){
450+ divide (nodes, xs, ys, zs, true , diag_balance);
451+ }
452+ // recurse into children
453+ for (int_t i = 0 ; i < (1 <<n_dim); ++i){
454+ children[i]->refine_image (nodes, image, shape_cells, xs, ys, zs, diag_balance);
455+ }
456+ }
457+ }
458+
421459void Cell::divide (node_map_t & nodes, double * xs, double * ys, double * zs, bool balance, bool diag_balance){
422460 // Gaurd against dividing a cell that is already at the max level
423461 if (level == max_level){
@@ -896,6 +934,18 @@ void Tree::refine_function(function test_func, bool diagonal_balance){
896934 roots[iz][iy][ix]->refine_func (nodes, test_func, xs, ys, zs, diagonal_balance);
897935};
898936
937+ void Tree::refine_image (double *image, bool diagonal_balance){
938+ int_t shape_cells[3 ];
939+ shape_cells[0 ] = nx/2 ;
940+ shape_cells[1 ] = ny/2 ;
941+ shape_cells[2 ] = nz/2 ;
942+ for (int_t iz=0 ; iz<nz_roots; ++iz)
943+ for (int_t iy=0 ; iy<ny_roots; ++iy)
944+ for (int_t ix=0 ; ix<nx_roots; ++ix)
945+ roots[iz][iy][ix]->refine_image (nodes, image, shape_cells, xs, ys, zs, diagonal_balance);
946+ }
947+
948+
899949void Tree::finalize_lists (){
900950 for (int_t iz=0 ; iz<nz_roots; ++iz)
901951 for (int_t iy=0 ; iy<ny_roots; ++iy)
0 commit comments