Skip to content

Commit

Permalink
#30 Order has been restored.
Browse files Browse the repository at this point in the history 
All functions of geometry now follow z,y,x, which corresponds to the rest of the code.
0700 now works on 770c_pag.
  • Loading branch information
@carljohnsen
carljohnsen committed Mar 5, 2024
1 parent fec3f47 commit 4711e5e
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 41 deletions.
50 changes: 25 additions & 25 deletions src/lib/cpp/cpu_seq/geometry.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ array<real_t, 3> center_of_mass(const input_ndarray<mask_type> &mask) {

print_timestamp("center_of_mass end");

return array<real_t, 3>{ rcmx, rcmy, rcmz };
return array<real_t, 3>{ rcmz, rcmy, rcmx };
}

void compute_front_mask(const input_ndarray<mask_type> solid_implant,
Expand Down Expand Up @@ -271,38 +271,38 @@ array<real_t,9> inertia_matrix(const input_ndarray<mask_type> &mask, const array
UNPACK_NUMPY(mask);

real_t
Ixx = 0, Ixy = 0, Ixz = 0,
Iyy = 0, Iyz = 0,
Izz = 0;
Izz = 0, Izy = 0, Izx = 0,
Iyy = 0, Iyx = 0,
Ixx = 0;

print_timestamp("inertia_matrix_serial start");

BLOCK_BEGIN(mask, reduction(+:Ixx, Iyy, Izz) reduction(-:Ixy,Ixz,Iyz)) {
BLOCK_BEGIN(mask, reduction(+:Izz, Iyy, Ixx) reduction(+:Izy,Izx,Iyx)) {

mask_type m = mask_buffer[flat_index];
real_t m = mask_buffer[flat_index];

// m guards this, and then branches are removed
//if (m != 0)
real_t
X = real_t(x) - cm[0],
Y = real_t(y) - cm[1],
Z = real_t(z) - cm[2];
Z = real(z) - cm[0],
Y = real(y) - cm[1],
X = real(x) - cm[2];

Ixx += m * (Y*Y + Z*Z);
Iyy += m * (X*X + Z*Z);
Izz += m * (X*X + Y*Y);
Ixy -= m * X*Y;
Ixz -= m * X*Z;
Iyz -= m * Y*Z;
Izz += m * (Y*Y + X*X);
Iyy += m * (Z*Z + X*X);
Ixx += m * (Z*Z + Y*Y);
Izy -= m * Z*Y;
Izx -= m * Z*X;
Iyx -= m * Y*X;

} BLOCK_END();

print_timestamp("inertia_matrix_serial end");

return array<real_t,9> {
Ixx, Ixy, Ixz,
Ixy, Iyy, Iyz,
Ixz, Iyz, Izz
Izz, Izy, Izx,
Izy, Iyy, Iyx,
Izx, Iyx, Ixx
};
}

Expand Down Expand Up @@ -376,21 +376,21 @@ void sample_plane(const input_ndarray<T> &voxels,

// X,Y,Z in micrometers; x,y,z in voxel index space
const real_t
X = cm[0] + u*u_axis[0] + v*v_axis[0],
Z = cm[0] + u*u_axis[0] + v*v_axis[0],
Y = cm[1] + u*u_axis[1] + v*v_axis[1],
Z = cm[2] + u*u_axis[2] + v*v_axis[2];
X = cm[2] + u*u_axis[2] + v*v_axis[2];

const real_t
z = Z / voxel_size,
x = X / voxel_size,
y = Y / voxel_size,
z = Z / voxel_size;
y = Y / voxel_size;

// printf("u,v = %g,%g -> %.1f,%.1f,%.1f -> %d, %d, %d\n",u,v,X,Y,Z,int(round(x)),int(round(y)),int(round(z)));

T value = 0;
std::array<float, 6> local_bbox = {0.5f, float(voxels_Nx)-0.5f, 0.5f, float(voxels_Ny)-0.5f, 0.5f, float(voxels_Nz)-0.5f};
if (in_bbox({{x,y,z}}, local_bbox))
value = (T) round(resample2x2x2<T>(voxels.data, {voxels_Nx, voxels_Ny, voxels_Nz}, {x, y, z}));
std::array<float, 6> local_bbox = {0.5f, float(voxels_Nz)-0.5f, 0.5f, float(voxels_Ny)-0.5f, 0.5f, float(voxels_Nx)-0.5f};
if (in_bbox({{z,y,x}}, local_bbox))
value = (T) round(resample2x2x2<T>(voxels.data, {voxels_Nz, voxels_Ny, voxels_Nx}, {z, y, x}));
// else
// fprintf(stderr,"Sampling outside image: x,y,z = %.1f,%.1f,%.1f, Nx,Ny,Nz = %ld,%ld,%ld\n",x,y,z,Nx,Ny,Nz);

Expand Down
2 changes: 1 addition & 1 deletion src/lib/cpp/include/geometry.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ float resample2x2x2(const T *voxels,
auto [Nz,Ny,Nx] = shape;

if (!in_bbox(index, {0.5f, float(Nx)-0.5f, 0.5f, float(Ny)-0.5f, 0.5f, float(Nz)-0.5f})) {
uint64_t voxel_index = uint64_t(floor(index[0]))*Nz*Ny + uint64_t(floor(index[1]))*Nx + uint64_t(floor(index[2]));
uint64_t voxel_index = uint64_t(floor(index[0]))*Ny*Nx + uint64_t(floor(index[1]))*Nx + uint64_t(floor(index[2]));
return voxels[voxel_index];
}

Expand Down
64 changes: 49 additions & 15 deletions src/processing_steps/0700_implant_FoR.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
sys.path.append(sys.path[0]+"/../")
from config.constants import *
from config.paths import hdf5_root, binary_root
from lib.cpp.cpu.geometry import center_of_mass, inertia_matrix, sample_plane
from lib.cpp.gpu.morphology import erode_3d_sphere as erode_3d, dilate_3d_sphere as dilate_3d
from lib.cpp.cpu_seq.geometry import center_of_mass, inertia_matrix, sample_plane
from lib.cpp.cpu.morphology import erode_3d_sphere as erode_3d, dilate_3d_sphere as dilate_3d
import matplotlib.pyplot as plt
from matplotlib.colors import colorConverter
import scipy as sp, scipy.ndimage as ndi, scipy.interpolate as interpolate, scipy.signal as signal
Expand Down Expand Up @@ -37,9 +37,9 @@ def close_3d(image, r):
(Nz,Ny,Nx) = image.shape
I1 = np.zeros((Nz+2*r,Ny+2*r,Nx+2*r), dtype=np.uint8)
I2 = np.zeros((Nz+2*r,Ny+2*r,Nx+2*r), dtype=np.uint8)
I1[r:-r,r:-r,r:-r] = image;
I1[r:-r,r:-r,r:-r] = image

dilate_3d(I1,r,I2);
dilate_3d(I1,r,I2)
erode_3d (I2,r,I1)

return I1[r:-r,r:-r,r:-r].astype(image.dtype)
Expand All @@ -48,9 +48,9 @@ def open_3d(image, r):
(Nz,Ny,Nx) = image.shape
I1 = np.zeros((Nz+2*r,Ny+2*r,Nx+2*r), dtype=np.uint8)
I2 = np.zeros((Nz+2*r,Ny+2*r,Nx+2*r), dtype=np.uint8)
I1[r:-r,r:-r,r:-r] = image;
I1[r:-r,r:-r,r:-r] = image

erode_3d (I1,r,I2);
erode_3d (I1,r,I2)
dilate_3d(I2,r,I1)

return I1[r:-r,r:-r,r:-r].astype(image.dtype)
Expand Down Expand Up @@ -119,7 +119,8 @@ def zyx_to_UVWp_transform():
daxis = {'z':np.array([-1,1,0]), 'y':np.array([0,0,1]), 'z2':np.array([-1.5,0,0])}

def figure_FoR_UVW(debug=True):
vol = vedo.Volume(implant, alpha=[0,0,0.05,0.2])
vol = vedo.Volume(implant)
vol.alpha([0,0,0.05,0.2])
u_arrow = vedo.Arrow(cm[::-1],cm[::-1]+1/np.sqrt(ls[0]/ls[2])*100*u_vec[::-1],c='r',s=0.7)
v_arrow = vedo.Arrow(cm[::-1],cm[::-1]+1/np.sqrt(ls[1]/ls[2])*100*v_vec[::-1],c='g',s=0.7)
w_arrow = vedo.Arrow(cm[::-1],cm[::-1]+100*w_vec[::-1],c='b',s=0.7)
Expand Down Expand Up @@ -235,7 +236,8 @@ def figure_FoR_cylinder(debug=True):
Vp_arrow = vedo.Arrow(Cp, UVW2xyz(cp+implant_radius*2*v_prime), c='g')
Wp_arrow = vedo.Arrow(Cp, UVW2xyz(cp+implant_radius*2*w_prime), c='b')

vol = vedo.Volume(implant,alpha=[0,0,0.05,0.1])
vol = vedo.Volume(implant)
vol.alpha([0,0,0.05,0.1])

pl = vedo.Plotter(offscreen=True, interactive=False,sharecam=False)
for axis in vaxis.keys():
Expand All @@ -251,7 +253,8 @@ def figure_FoR_cylinder(debug=True):
vedo.show([vol,cylinder,Up_arrow,Vp_arrow,Wp_arrow],interactive=True)

def figure_FoR_voxels(name,voxels,debug=True):
vol = vedo.Volume(voxels,alpha=[0,0,0.05,0.1])
vol = vedo.Volume(voxels)
vol.alpha([0,0,0.05,0.1])

pl = vedo.Plotter(offscreen=True, interactive=False,sharecam=False)
for axis in vaxis.keys():
Expand Down Expand Up @@ -289,11 +292,11 @@ def figure_FoR_voxels(name,voxels,debug=True):
voxels = np.fromfile(f"{binary_root}/voxels/{scale}x/{sample}.uint16",dtype=np.uint16).reshape(implant.shape)

plt.imshow(implant[implant.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy.png')
plt.imshow(implant[:,implant.shape[0]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz.png')
plt.imshow(implant[:,:,implant.shape[0]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz.png')
plt.imshow(implant[:,implant.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz.png')
plt.imshow(implant[:,:,implant.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz.png')
plt.imshow(voxels[voxels.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/voxels-sanity-xy.png')
plt.imshow(voxels[:,voxels.shape[0]//2,:]); plt.savefig(f'{image_output_dir}/voxels-sanity-xz.png')
plt.imshow(voxels[:,:,voxels.shape[0]//2]); plt.savefig(f'{image_output_dir}/voxels-sanity-yz.png')
plt.imshow(voxels[:,voxels.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/voxels-sanity-xz.png')
plt.imshow(voxels[:,:,voxels.shape[2]//2]); plt.savefig(f'{image_output_dir}/voxels-sanity-yz.png')

nz,ny,nx = implant.shape

Expand Down Expand Up @@ -456,7 +459,10 @@ def UVW2xyz(p):
figure_FoR_voxels("solid_implant",solid_implant,verbose >= 2)

back_mask = (Ws<0)
front_mask = largest_cc_of((Ws>50)*(~solid_implant))#*(thetas>=theta_from)*(thetas<=theta_to)
front_mask = largest_cc_of((Ws>50)&(~solid_implant))#*(thetas>=theta_from)*(thetas<=theta_to)
plt.imshow(front_mask[front_mask.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-front_mask.png')
plt.imshow(front_mask[:,front_mask.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-front_mask.png')
plt.imshow(front_mask[:,:,front_mask.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-front_mask.png')

# back_part = voxels*back_mask

Expand Down Expand Up @@ -516,24 +522,36 @@ def UVW2xyz(p):
group_name="implant_solid",
datasets={"mask":solid_implant},
attributes={"sample":sample,"scale":scale,"voxel_size":voxel_size})
plt.imshow(solid_implant[solid_implant.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-solid.png')
plt.imshow(solid_implant[:,solid_implant.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-solid.png')
plt.imshow(solid_implant[:,:,solid_implant.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-solid.png')

if verbose >= 1: print(f"Saving implant_shell mask to {output_dir}/{sample}.h5")
update_hdf5_mask(f"{output_dir}/{sample}.h5",
group_name="implant_shell",
datasets={"mask":implant_shell_mask},
attributes={"sample":sample,"scale":scale,"voxel_size":voxel_size})
plt.imshow(implant_shell_mask[implant_shell_mask.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-shell.png')
plt.imshow(implant_shell_mask[:,implant_shell_mask.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-shell.png')
plt.imshow(implant_shell_mask[:,:,implant_shell_mask.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-shell.png')

if verbose >= 1: print(f"Saving cut_cylinder_air mask to {output_dir}/{sample}.h5")
update_hdf5_mask(f"{output_dir}/{sample}.h5",
group_name="cut_cylinder_air",
datasets={"mask":back_mask},
attributes={"sample":sample,"scale":scale,"voxel_size":voxel_size})
plt.imshow(back_mask[back_mask.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-back.png')
plt.imshow(back_mask[:,back_mask.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-back.png')
plt.imshow(back_mask[:,:,back_mask.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-back.png')

if verbose >= 1: print(f"Saving cut_cylinder_bone mask to {output_dir}/{sample}.h5")
update_hdf5_mask(f"{output_dir}/{sample}.h5",
group_name="cut_cylinder_bone",
datasets={"mask":front_mask},
attributes={"sample":sample, "scale":scale, "voxel_size":voxel_size})
plt.imshow(front_mask[front_mask.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-front.png')
plt.imshow(front_mask[:,front_mask.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-front.png')
plt.imshow(front_mask[:,:,front_mask.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-front.png')

if verbose >= 1: print(f"Computing bone region")
hist, bins = np.histogram(front_part, 256)
Expand All @@ -546,6 +564,10 @@ def UVW2xyz(p):
if verbose >= 1: print(f"p1, p2 = ({p1,bins[p1]}), ({p2,bins[p2]}); midpoint = {midpoint}")

bone_mask1 = front_part > midpoint
plt.imshow(bone_mask1[bone_mask1.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-bone1.png')
plt.imshow(bone_mask1[:,bone_mask1.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-bone1.png')
plt.imshow(bone_mask1[:,:,bone_mask1.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-bone1.png')

closing_diameter, opening_diameter = 400, 300 # micrometers
closing_voxels = 2*int(round(closing_diameter/(2*voxel_size))) + 1 # Scale & ensure odd length
opening_voxels = 2*int(round(opening_diameter/(2*voxel_size))) + 1 # Scale & ensure odd length
Expand All @@ -554,10 +576,19 @@ def UVW2xyz(p):
bone_region_mask = close_3d(bone_mask1, closing_voxels//2)

for i in tqdm.tqdm(range(1),f"Opening with sphere of diameter {opening_diameter} micrometers, {opening_voxels} voxels.\n"):
bone_region_mask &= ~solid_implant #~open_3d(implant_shell_mask, opening_voxels)
bone_region_mask &= (~solid_implant).astype(bool) #~open_3d(implant_shell_mask, opening_voxels)
bone_region_mask = open_3d(bone_region_mask,opening_voxels//2)

bone_region_mask = largest_cc_of(bone_region_mask)
plt.imshow(bone_region_mask[bone_region_mask.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-bone.png')
plt.imshow(bone_region_mask[:,bone_region_mask.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-bone.png')
plt.imshow(bone_region_mask[:,:,bone_region_mask.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-bone.png')

if verbose >= 1: print(f"Saving bone_region mask to {output_dir}/{sample}.h5")
update_hdf5_mask(f"{output_dir}/{sample}.h5",
group_name="bone_region",
datasets={"mask":bone_region_mask},
attributes={"sample":sample, "scale":scale, "voxel_size":voxel_size})
except:
if verbose >= 1: print(f"Wasnt able to separate into resin and bone region. Assuming all is bone region.")
bone_region_mask = front_mask
Expand All @@ -567,3 +598,6 @@ def UVW2xyz(p):
group_name="bone_region",
datasets={"mask":bone_region_mask},
attributes={"sample":sample, "scale":scale, "voxel_size":voxel_size})
plt.imshow(bone_region_mask[bone_region_mask.shape[0]//2,:,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xy-bone.png')
plt.imshow(bone_region_mask[:,bone_region_mask.shape[1]//2,:]); plt.savefig(f'{image_output_dir}/implant-sanity-xz-bone.png')
plt.imshow(bone_region_mask[:,:,bone_region_mask.shape[2]//2]); plt.savefig(f'{image_output_dir}/implant-sanity-yz-bone.png')

2 comments on commit 4711e5e

@carljohnsen
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should be #38, not #30

@jamesavery
Copy link
Owner

@jamesavery jamesavery commented on 4711e5e Mar 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Grossartig! Ordnung muss sein!

Please sign in to comment.