From 5d9ae4db3a06d9c12803fefd3e34b78244c8bed0 Mon Sep 17 00:00:00 2001
From: Jean-Yves TINEVEZ <jean-yves.tinevez@pasteur.fr>
Date: Thu, 16 Jul 2020 11:46:31 +0200
Subject: [PATCH] Externalize curvature calculation and put the right units.
---
src/@deproj/compute_curvatures.m | 48 ++++++++++++++++++++++++++++++++
src/@deproj/deproj.m | 4 +++
src/@deproj/from_heightmap.m | 28 ++++---------------
src/@deproj/plot_curvatures.m | 3 +-
4 files changed, 60 insertions(+), 23 deletions(-)
create mode 100644 src/@deproj/compute_curvatures.m
diff --git a/src/@deproj/compute_curvatures.m b/src/@deproj/compute_curvatures.m
new file mode 100644
index 0000000..e02e83c
--- /dev/null
+++ b/src/@deproj/compute_curvatures.m
@@ -0,0 +1,48 @@
+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
+
diff --git a/src/@deproj/deproj.m b/src/@deproj/deproj.m
index d8cd05e..5b02320 100644
--- a/src/@deproj/deproj.m
+++ b/src/@deproj/deproj.m
@@ -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.
diff --git a/src/@deproj/from_heightmap.m b/src/@deproj/from_heightmap.m
index b20fc6a..9659342 100644
--- a/src/@deproj/from_heightmap.m
+++ b/src/@deproj/from_heightmap.m
@@ -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.
diff --git a/src/@deproj/plot_curvatures.m b/src/@deproj/plot_curvatures.m
index 051242d..46fc54e 100644
--- a/src/@deproj/plot_curvatures.m
+++ b/src/@deproj/plot_curvatures.m
@@ -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 );
--
GitLab