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

Externalize curvature calculation and put the right units.

parent d0cd038e
function [ curvMean, curvGauss, curvK1, curvK2 ] = compute_curvatures( H, object_scale, pixel_size, voxel_depth, invert_z )
%COMPUTE_CURVATURES Compute local curvature from the smoothed height-map.
% Ref: http://users.vcnet.com/simonp/curvature.pdf
if nargin < 5
invert_z = false;
end
if nargin < 4
voxel_depth = 1;
end
if nargin < 3
pixel_size = 1;
end
if nargin >= 2 && ~isnan( object_scale) && ~isempty( object_scale )
% We need to smooth the height-map over the scale of several cells.
Hs = imgaussfilt( H, 3 * object_scale );
else
Hs = H;
end
% Physical units.
Hs = Hs * voxel_depth;
if invert_z
Hs = -Hs;
end
[ Hx , Hy ] = gradient( Hs, pixel_size );
[ Hxx, Hxy ] = gradient( Hx, pixel_size );
[ ~ , Hyy ] = gradient( Hy, pixel_size );
% Gaussian curvature.
Nk = ( 1. + Hx .^ 2 + Hy .^ 2 );
curvGauss = ( Hxx .* Hyy - Hxy .^ 2 ) ./ ( Nk .^ 2 );
% Mean curvature.
Dk1 = ( 1. + Hx .^ 2 ) .* Hyy;
Dk2 = ( 1. + Hy .^ 2 ) .* Hxx;
Dk3 = - 2. * Hx .* Hy .* Hxy;
curvMean = ( Dk1 + Dk2 + Dk3 ) ./ ( 2 * Nk .^ 1.5 );
% Principle curvatures.
curvK1 = curvMean + sqrt( curvMean .* curvMean - curvGauss );
curvK2 = curvMean - sqrt( curvMean .* curvMean - curvGauss );
end
......@@ -87,6 +87,10 @@ classdef deproj
% Returns the seismic colormap.
cmap = cmap_seismic();
% Compute local curvature from the smoothed height-map.
[ curvMean, curvGauss, curvK1, curvK2 ] = compute_curvatures( H, object_scale, pixel_size, voxel_depth, invert_z )
end
%% Private static methods: utilities.
......
......@@ -139,28 +139,12 @@ function obj = from_heightmap( ...
% Ref: http://users.vcnet.com/simonp/curvature.pdf
fprintf('Compuing tissue local curvature.\n' )
% We need to smooth the height-map over the scale of several cells.
Hs2 = imgaussfilt( Hs1, 3 * object_scale );
if invert_z
Hs2 = -Hs2;
end
[ Hx , Hy ] = gradient( Hs2 );
[ Hxx, Hxy ] = gradient( Hx );
[ ~ , Hyy ] = gradient( Hy );
% Gaussian curvature.
Nk = ( 1. + Hx .^ 2 + Hy .^ 2 );
curvGauss = ( Hxx .* Hyy - Hxy .^ 2 ) ./ ( Nk .^ 2 );
% Mean curvature.
Dk1 = ( 1. + Hx .^ 2 ) .* Hyy;
Dk2 = ( 1. + Hy .^ 2 ) .* Hxx;
Dk3 = - 2. * Hx .* Hy .* Hxy;
curvMean = ( Dk1 + Dk2 + Dk3 ) ./ ( 2 * Nk .^ 1.5 );
% Principle curvatures.
curvK1 = curvMean + sqrt( curvMean .* curvMean - curvGauss );
curvK2 = curvMean - sqrt( curvMean .* curvMean - curvGauss );
[ curvMean, curvGauss, curvK1, curvK2 ] = deproj.compute_curvatures( ...
Hs1, ...
object_scale, ...
pixel_size, ...
voxel_depth, ...
invert_z );
%% Create epicell instances.
......
......@@ -43,7 +43,8 @@ function [ hf, ax1, ax2, ax3 ] = plot_curvatures( obj, scale_bar_length )
colormap( ax1, cmap )
colormap( ax2, cmap )
colormap( ax3, cmap )
colorbar(ax3, 'Location', 'EastOutside' )
c = colorbar(ax3, 'Location', 'EastOutside' );
c.Label.String = sprintf('Curvature (1/%s)', obj.units );
add_plot_scalebar( obj, scale_bar_length, ax3 );
......
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