Hi Yan,
great code. I have one suggestion though: In your CV reinit process, you implemented a very efficient distance transform. However, for performance reasons, you don't compute the values for the edges of the matrix (always 0). Thus I would suggest, you pad the input matrix with 0s and then un-pad the output matrix of this step:
u0 = y_binary_boundary_detection(uint8(u>0));
u0 = padarray(u0, [1 1 1]); % padding
u0 = ac_distance_transform_3d(u0);
u0 = u0(2:end-1, 2:end-1, 2:end-1); % un-padding
u = u0.*sign(u);
Cheers, Christian