Commit d6dc1e60 authored by Jean-Yves TINEVEZ's avatar Jean-Yves TINEVEZ
Browse files

Refine the height-map processing.

parent 326773c8
......@@ -117,15 +117,18 @@ classdef deproj
% Create a height-map from a mesh.
[ H, min_y, min_x ] = mesh_to_heightmap( V, pixel_size )
% Prepare the height-map image according to user settings.
H = prepare_heightmap( H, voxel_depth, smooth_scale, invert_z, inpaint_zeros, prune_zeros )
% Returns the Z position of points taken from a processed height-map.
z_coords = get_z( P, H, pixel_size )
end
%% Private static methods: utilities.
methods ( Access = private, Hidden = true, Static = true )
% Returns the Z position of points taken from a height-map.
z_coords = get_z( P, H, pixel_size, voxel_depth )
% Returns the cells from a BW image with ridges.
[ objects, junction_graph ] = mask_to_objects( I, downsample )
......
......@@ -50,10 +50,10 @@ function obj = from_heightmap( ...
med_area = median(arrayfun( @(s) polyarea( s.boundary(:,1), s.boundary(:,2) ), ...
objects ) ); % pixels^2
object_scale = 2 * sqrt( med_area )/ pi; % pixel units
smooth_scale = 2 * sqrt( med_area )/ pi; % pixel units
fprintf('Typical object scale: %.1f pixels or %.2f µm.\n', ...
object_scale, object_scale * pixel_size )
smooth_scale, smooth_scale * pixel_size )
%% Scale to physical coordinates.
......@@ -71,25 +71,16 @@ function obj = from_heightmap( ...
fprintf('Collecting Z coordinates.\n' )
tic
if inpaint_zeros
mask = H == 0;
H = regionfill( H, mask );
end
if prune_zeros
H( H == 0 ) = NaN;
end
% Smooth the height-map over a scale smaller than a cell.
Hs1 = imgaussfilt( H, object_scale );
H = deproj.prepare_heightmap( H, voxel_depth, smooth_scale, invert_z, inpaint_zeros, prune_zeros );
% For junction.
z_junction = deproj.get_z( junction_graph.Nodes.Centroid, Hs1, pixel_size, voxel_depth );
z_junction = deproj.get_z( junction_graph.Nodes.Centroid, H, pixel_size );
junction_graph.Nodes.Centroid = [ junction_graph.Nodes.Centroid z_junction ];
% For objects.
for i = 1 : n_objects
z_obj = deproj.get_z( objects( i ).boundary, Hs1, pixel_size, voxel_depth );
z_obj = deproj.get_z( objects( i ).boundary, H, pixel_size );
objects( i ).boundary = [ objects( i ).boundary z_obj ];
objects( i ).center(3) = mean( z_obj );
end
......@@ -120,31 +111,17 @@ function obj = from_heightmap( ...
n_objects = numel( objects );
end
%% Invert Z position for plotting.
if invert_z
max_z = max( z_junction );
z_junction = max_z - z_junction;
junction_graph.Nodes.Centroid = [ junction_graph.Nodes.Centroid z_junction ];
for i = 1 : n_objects
objects( i ).boundary( :, 3 ) = max_z - objects( i ).boundary( :, 3 );
objects( i ).center( 3 ) = max_z - objects( i ).center( 3 );
end
end
end
%% Compute local curvature from the smoothed height-map.
% Ref: http://users.vcnet.com/simonp/curvature.pdf
fprintf('Computing tissue local curvature.\n' )
[ curvMean, curvGauss, curvK1, curvK2 ] = deproj.compute_curvatures( ...
Hs1, ...
object_scale, ...
H, ...
smooth_scale, ...
pixel_size, ...
voxel_depth, ...
invert_z );
voxel_depth );
%% Create epicell instances.
......
function z_coords = get_z( P, H, pixel_size, voxel_depth )
function z_coords = get_z( P, H, pixel_size )
%GET_Z Returns the Z position of points taken from a height-map.
% - P is a Nx2 list of points, in physical coordinates.
% - H is the height map, encoding the z plane of interest for all X & Y.
% - pixel_size: convert pixel coordinates to physical coordinates.
% - voxel_size: convert plane of interest to physical coordinates.
% Returns the Z coordinates vector in physical coordinates.
% Returns the Z coordinates vector in physical coordinates.
%% Yield coordinates.
Pp = round( P / pixel_size ); % pixel coordinates.
......@@ -19,7 +20,7 @@ function z_coords = get_z( P, H, pixel_size, voxel_depth )
yp( yp > height ) = height;
xy_ind = sub2ind( size(H), yp, xp );
z_coords = H( xy_ind ) * voxel_depth;
z_coords = H( xy_ind );
end
function H = prepare_heightmap( H, voxel_depth, smooth_scale, invert_z, inpaint_zeros, prune_zeros )
%PREPARE_HEIGHTMAP Prepare the height-map image according to user settings.
if nargin < 6
prune_zeros = true;
end
if nargin < 5
inpaint_zeros = true;
end
if nargin < 4
invert_z = false;
end
if nargin < 3
smooth_scale = -1.;
end
if nargin < 2
voxel_depth = 1.;
end
%% Process height-map.
if inpaint_zeros
mask = H == 0;
H = regionfill( H, mask );
end
if prune_zeros
H( H == 0 ) = NaN;
end
if smooth_scale > 0
H = imgaussfilt( H, smooth_scale );
end
if invert_z
H = max(H(:)) - H;
end
H = H * voxel_depth;
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment