From d94757dc96cb04697364608148e0f8361f1bc3d8 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 9 Apr 2021 13:47:11 +0100 Subject: [PATCH 01/30] working on setting up interpolators --- createInterpolator.m | 46 ++++++++++++++++ getConductionVelocity.m | 114 ++++++++++++++++++++++++++++++++------- openEpDataInterpolator.m | 83 ++++++++++++++++++++++++++++ 3 files changed, 224 insertions(+), 19 deletions(-) create mode 100644 createInterpolator.m create mode 100644 openEpDataInterpolator.m diff --git a/createInterpolator.m b/createInterpolator.m new file mode 100644 index 0000000..dd99f8b --- /dev/null +++ b/createInterpolator.m @@ -0,0 +1,46 @@ +function int = createInterpolator( type, options ) +% CREATEINTERPOLATOR Creates an interpolation scheme for use with OpenEP +% data +% +% Usage: +% int = createInterpolator( type ) +% +% Where: +% type - the type of interpolator required which may be one of, +% `'scatteredInterpolant'`, `'radialbasis'` or `'localsmoothing'` +% int - a function handle to the interpolator +% options - a structure with the following fields: +% options.method +% options.extrapolationmethod +% options.distancethreshold +% +% CREATEINTERPOLATOR creates and interpolation scheme for use with OpenEP +% data +% +% Author: Steven Williams (2021) (Copyright) +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% cvdata = getConductionVelocity( userdata ); +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +switch(type) + + case 'scatteredInterpolant' + + method = options.method; + extrapolationMethod = options.extrapolationMethod; + distanceThreshold = options.distanceThreshold; + + case 'radialbasis' + + case 'localsmoothing' + +end \ No newline at end of file diff --git a/getConductionVelocity.m b/getConductionVelocity.m index ab94152..cfd63ab 100644 --- a/getConductionVelocity.m +++ b/getConductionVelocity.m @@ -1,22 +1,52 @@ -function cvdata = getConductionVelocity( userdata ) -% GETCONDUCTIONVELOCITY Returns the conduction velocity map of the -% chamber +function [cv, cvX, interpCv] = getConductionVelocity( userdata, varargin ) +% GETCONDUCTIONVELOCITY Returns the conduction velocity map of the chamber % % Usage: % cvdata = getConductionVelocity( userdata ) +% cvdata = getConductionVelocity( userrdata, interpolator) +% cvdata = getConductionVelocity( ... , param, value ) % Where: -% userdata - see importcarto_mem -% cvdata - the conduction velocities, in m/s -% -% GETCONDUCTIONVELOCITY Calculate conduction velocities by calculating -% gradients of interpolated local activation times. GETCONDUCTIONVELOCITY -% makes use of a modified version of "Scattered Data Interpolation and -% Approximation using Radial Base Functions" available from the Matlab -% FileExchange: Alex Chirokov (2020). Scattered Data Interpolation and -% Approximation using Radial Base Functions -% (https://www.mathworks.com/matlabcentral/fileexchange/10056-scattered-data-interpolation-and-approximation-using-radial-base-functions), MATLAB Central File Exchange. Retrieved November 24, 2020. -% -% Author: Steven Williams (2020) (Copyright) +% userdata - see importcarto_mem +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% GETCONDUCTIONVELOCITY accepts the following parameter-value pairs +% `'method'` `'triangulation'` | `'cosinefit'` | `{'radialbasis'}` | +% `'omnipole'` | `'gradient'` | `'eikonal'` +% `'interpolator'` `{'scatteredInterpolant'}` | `'radialbasis'` | +% `'localsmoothing'` + +% GETCONDUCTIONVELOCITY Returns the conduction velocity data of the chamber. +% Five methods for calculating conduction velocity are provided as +% described below. The conduction velocity data interpolated over the +% surface of the shell is also returned. The desired interpolation scheme +% can be set by passing in an OpenEPInterpolator as the second argument. +% +% Conduction velocity methods +% --------------------------- +% (*) triangulation +% description goes here +% +% (*) cosinefit +% description goes here +% +% (*) radialbasis +% description goes here +% +% (*) omnipole +% description goes here +% +% (*) gradient +% description goes here +% +% (*) eikonal +% description goes here +% +% Author: Steven Williams (2021) (Copyright) % SPDX-License-Identifier: Apache-2.0 % % Modifications - @@ -30,9 +60,55 @@ % code % --------------------------------------------------------------- -lats = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); -X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); -[~,~,~,~,cvdata] = RBFConductionVelocity(lats', X', userdata.surface.triRep.X'); -cvdata = cvdata'; +nStandardArgs = 1; % UPDATE VALUE +method = 'radialbasis'; +interpolator = 'scatteredinterpolant'; +if nargin > nStandardArgs + if ischar(varargin{2}) + for i = 1:2:nargin-nStandardArgs + switch varargin{i} + case 'method' + method = varargin{i+1}; + case 'interpolator' + interpolator = varargin{i+1}; + end + end + else + interpolator = varargin{2}; + end +end +% first create an interpolator +if ischar(interpolator) + int = createInterpolator(interpolator); +else + int = interpolator; +end + + switch lower(method) + case 'triangulation' + [cv, cvX, interpCv] = doCvMapping_Triangulation(userdata, int); + + case 'cosinefit' + [cv, cvX, interpCv] = doCvMapping_CosineFit(userdata); + + case 'radialbasis' + [cv, cvX, interpCv] = doCvMapping_RadialBasis(userdata); + + % lats = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); + % X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); + % [~,~,~,~,cvdata] = RBFConductionVelocity(lats', X', userdata.surface.triRep.X'); + % cvdata = cvdata'; + + case 'omnipole' + [cv, cvX, interpCv] = doCvMapping_Omnipole(userdata, int); + + case 'gradient' + [cv, cvX, interpCv] = doCvMapping_Gradient(userdata, int); + + case 'eikonal' + [cv, cvX, interpCv] = doCvMapping_Eikonal(userdata, int); + + end + end \ No newline at end of file diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m new file mode 100644 index 0000000..64ad2d9 --- /dev/null +++ b/openEpDataInterpolator.m @@ -0,0 +1,83 @@ +classdef openEpDataInterpolator + % OPENEPDATAINTERPOLATOR Performs data interpolation for OpenEP + + properties + method + interMethod + exterMethod + distanceThreshold + smoothingLength + end + + methods + % Contructor + function int = openEpDataInterpolator(varargin) + + % specify default values + int.interMethod = 'linear'; + int.exterMethod = 'linear'; + int.distanceThreshold = Inf; + int.smoothingLength = 10; + int.method = 'scatteredinterpolant'; + + if nargin==1 + int.method = varargin{1}; + end + if nargin==2 + options = varargin{2}; + if isfield(options, 'interMethod') + if ~isempty(options.interMethod) + int.interMethod = options.interMethod; + end + end + if isfield(options, 'exterMethod') + if ~isempty(options.exterMethod) + int.exterMethod = options.exterMethod; + end + end + if isfield(options, 'distanceThreshold') + if ~isempty(options.distanceThreshold) + int.distanceThreshold = options.distanceThreshold; + end + end + end + end + + % Interpolation + function d1 = interpolate(int, x0,d0,x1) + + % remove any data with NaNs + tempData = d0; + tempCoords = x0; + iNaN = isnan(tempData); + tempData(iNaN) = []; + tempCoords(iNaN,:) = []; + d0 = tempData; + x0 = tempCoords; + + % perform interpolation + switch lower(int.method) + case 'scatteredinterpolant' + if numel(d0) < 4 + error('OPENEP/OPENEPDATAINTERPOLATOR: Cannot perform interpolation with less than four data points'); + end + F = scatteredInterpolant(x0(:,1), x0(:,2), x0(:,3) ... + , d0, int.interMethod, int.exterMethod); + d1 = F(x1); + + case 'localsmoothing' + + case 'radialbasis' + end + + % remove interpolated data that is too far from real data + id = knnsearch(x0, x1); + cPts = x0(id,:); %c for closest + d = distBetweenPoints(cPts, x1); + thresholdDistance = zeros(size(d)); + thresholdDistance(d>int.distanceThreshold) = 1; + nanset = logical(thresholdDistance); + d1(nanset) = NaN; + end + end +end \ No newline at end of file From 85182e9f7d3ef504111076fe2ed266757a471b1b Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 9 Apr 2021 14:20:24 +0100 Subject: [PATCH 02/30] Finished creating an initial version of openEPDataInterpolator.m --- openEpDataInterpolator.m | 50 ++++++++++++++++++++++++++--- private/localSmoothing.m | 68 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 114 insertions(+), 4 deletions(-) create mode 100644 private/localSmoothing.m diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m index 64ad2d9..3aa8c8f 100644 --- a/openEpDataInterpolator.m +++ b/openEpDataInterpolator.m @@ -3,10 +3,14 @@ properties method + end + + properties (Access = private) interMethod exterMethod distanceThreshold smoothingLength + fillWith end methods @@ -18,6 +22,8 @@ int.exterMethod = 'linear'; int.distanceThreshold = Inf; int.smoothingLength = 10; + int.fillWith = NaN; + int.method = 'scatteredinterpolant'; if nargin==1 @@ -27,12 +33,18 @@ options = varargin{2}; if isfield(options, 'interMethod') if ~isempty(options.interMethod) - int.interMethod = options.interMethod; + if ~strcmpi(int.method, 'scatteredInterpolant') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting interMethod property has no effect for interpolation method ' int.method]); + int.interMethod = options.interMethod; + end end end if isfield(options, 'exterMethod') if ~isempty(options.exterMethod) - int.exterMethod = options.exterMethod; + if ~strcmpi(int.method, 'scatteredInterpolant') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting exterMethod property has no effect for interpolation method ' int.method]); + int.exterMethod = options.exterMethod; + end end end if isfield(options, 'distanceThreshold') @@ -40,6 +52,22 @@ int.distanceThreshold = options.distanceThreshold; end end + if isfield(options, 'smoothingLength') + if ~isempty(options.smoothingLength) + if ~strcmpi(int.method, 'localSmoothing') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) + int.smoothingLength = options.smoothingLength; + end + end + end + if isfield(options, 'fillWith') + if ~isempty(options.fillWith) + if ~strcmpi(int.method, 'localSmoothing') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) + int.fillWith = options.fillWith; + end + end + end end end @@ -58,16 +86,30 @@ % perform interpolation switch lower(int.method) case 'scatteredinterpolant' + + % error checking if numel(d0) < 4 error('OPENEP/OPENEPDATAINTERPOLATOR: Cannot perform interpolation with less than four data points'); end - F = scatteredInterpolant(x0(:,1), x0(:,2), x0(:,3) ... - , d0, int.interMethod, int.exterMethod); + + % interpolation + F = scatteredInterpolant(x0(:,1), x0(:,2), x0(:,3), d0, int.interMethod, int.exterMethod); d1 = F(x1); case 'localsmoothing' + % error checking + % % % AC: is there any appropriate error checking to do + % here? + + % interpolation + d1 = localSmoothing(x0, d0, x1, int.smoothingLength, int.fillWith); + case 'radialbasis' + + % error checking + + % interpolation end % remove interpolated data that is too far from real data diff --git a/private/localSmoothing.m b/private/localSmoothing.m new file mode 100644 index 0000000..69743ec --- /dev/null +++ b/private/localSmoothing.m @@ -0,0 +1,68 @@ +function [f_dash, df_dash] = localSmoothing(x, f, x_dash, smoothingLength, fillWith) +% ------------ +% AC 8/2/21: local smoothing or piecewise interpolation, and locally +% centered moving least-squares gradient +% ------------ +% x is the points at which field values f are defined. +% size(x) = [m, 3], where m is the number of points. +% f are the values at the points x. size(f) = [m, 1]. +% x_dash are the points at which to query f. size(x_dash) = [n, 3]. +% smoothingLength is a scalar smoothing length. +% size(smoothingLength) = [1, 1] +% f_dash is the estimated function at points x_dash. size(f_dash) = +% [n, 1] +% df_dash is the estimted function gradient at x_dash. size(df_dash) = [n, +% 3] +% ------------ + +f_dash = zeros(length(x_dash), 1); +df_dash = zeros(length(x_dash), 3); +[Idx, dists] = rangesearch(x, x_dash, smoothingLength); +no_smooth = cellfun(@isempty, Idx); +if strcmp(fillWith, 'nearestneighbour') + nnIdx = knnsearch(x, x_dash(no_smooth, :)); + f_dash(no_smooth) = f(nnIdx); +else + f_dash(no_smooth) = NaN; +end + +kern = @(x, h) exp(-(x./h).^2); +% kern = @cubicSpline; + +for i=1:numel(Idx) + + if ~isempty(Idx{i}) + + dx = x(Idx{i},:) - x_dash(i, :); + + wi = kern(dists{i}', smoothingLength); + psi = wi./sum(wi); + + Xi = zeros(3); + for j = 1:size(dx, 1) + Xi = Xi + psi(j)*dx(j, :)'*dx(j, :); + end + bi = pinv(Xi); + + f_dash(i) = sum(f(Idx{i}).*psi); + + for j=1:size(dx, 1) + df_dash(i, :) = df_dash(i, :) + (f(Idx{i}(j)) - f_dash(i))*dx(j,:)*bi*psi(j); + end + + end + +end + + + function w = cubicSpline(d, smoothingLength) + q = d./smoothingLength; + d1 = 1 - 1.5*q.^2.*(1 - q./2); + d2 = 0.25*(2 - q).^3; + w = zeros(size(d)); + w(q<=1) = d1(q<=1); + w(q>1 & q<=2) = d2(q>1 & q<=2); + end + + +end \ No newline at end of file From 8aca9cbef73fdac72b2c61528f71dc327cccf5c1 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 9 Apr 2021 16:00:54 +0100 Subject: [PATCH 03/30] Added getVoltages, modified getVertices and finished setting up the interpolator class for scattered and local smoothing --- getVertices.m | 26 ++++++++++++-- getVoltages.m | 49 ++++++++++++++++++++++++++ openEpDataInterpolator.m | 74 +++++++++++++++++++++++++++++++++++----- 3 files changed, 138 insertions(+), 11 deletions(-) create mode 100644 getVoltages.m diff --git a/getVertices.m b/getVertices.m index f324acb..7f1fced 100644 --- a/getVertices.m +++ b/getVertices.m @@ -1,4 +1,4 @@ -function [vertices, isVertUsed] = getVertices( userdata ) +function [vertices, isVertUsed] = getVertices( userdata, varargin ) % GETVERTICES Returns the vertices referenced by userdata % % Usage: @@ -8,12 +8,18 @@ % vertices - all the vertices % isVertUsed - whether the vertex is referenced by the triangulation % -% GETVERTICES Returns the vertices referenced by userdata +% GETVERTICES accepts the following parameter-value pairs +% 'used' {false} | true +% +% GETVERTICES Returns the vertices referenced by userdata. If the property +% `used` is set to `true` then only vertices in the triangulation that are +% used by the triangulation are returned. % % Author: Steven Williams (2020) (Copyright) % SPDX-License-Identifier: Apache-2.0 % % Modifications - +% SW 09-04-2021: add routine to only return used vertices % % Info on Code Testing: % --------------------------------------------------------------- @@ -24,7 +30,21 @@ % code % --------------------------------------------------------------- +nStandardArgs = 1; % UPDATE VALUE +used = false; +if nargin > nStandardArgs + for i = 1:2:nargin-nStandardArgs + switch varargin{i} + case 'used' + used = varargin{i+1}; + end + end +end + vertices = userdata.surface.triRep.X; [~, isVertUsed] = repack(getMesh(userdata)); -end +if used + vertices = vertices(isVertUsed,:); + isVertUsed = true(size(vertices)); +end \ No newline at end of file diff --git a/getVoltages.m b/getVoltages.m new file mode 100644 index 0000000..b2db221 --- /dev/null +++ b/getVoltages.m @@ -0,0 +1,49 @@ +function [voltages] = getVoltages( userdata, varargin ) +% GETVOLTAGES Returns the voltages associated with each mapping point +% +% Usage: +% [voltages] = getVoltages( userdata ) +% Where: +% userdata - see importcarto_mem +% voltages - all the voltages +% +% GETVOLTAGES accepts the following parameter-value pairs +% 'type' {'bip'} | 'uni' +% +% GETVOLTAGES Returns the voltage values associated with each mapping +% point. Bipolar or unipolar voltages can be specified by setting the +% `type` property. +% +% Author: Steven Williams (2021) (Copyright) +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% voltages = getVoltages( userdata, 'type', 'bip' ) +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +nStandardArgs = 1; % UPDATE VALUE +type = 'bip'; +if nargin > nStandardArgs + for i = 1:2:nargin-nStandardArgs + switch varargin{i} + case 'type' + type = varargin{i+1}; + end + end +end + +switch lower(type) + case 'bip' + voltages = userdata.electric.voltages.bipolar; + case 'uni' + voltages = userdata.electric.voltages.unipolar; +end + +end \ No newline at end of file diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m index 3aa8c8f..2bc8a00 100644 --- a/openEpDataInterpolator.m +++ b/openEpDataInterpolator.m @@ -1,11 +1,69 @@ classdef openEpDataInterpolator - % OPENEPDATAINTERPOLATOR Performs data interpolation for OpenEP + % OPENEPDATAINTERPOLATOR Creates objects for performing spatial + % interpolation for OpenEP data + % + % Usage (Constructor): + % int = openEpDataInterpolator() + % int = openEpDataInterpolator(method) + % int = openEpDataInterpolator(method, options) + % + % Where: + % method - the method to use for interpolation, can be + % 'scatteredInterpolant', 'localSmoothing' or 'radialBasis' + % options - a structure containing options pertaining to each + % interpolation method. Options which are not relevant to a + % particular interpolation method will be ignored. + % .interMethod + % - The interpolation method to use with + % scatteredInterpolant. See `help scatteredInterpolant` + % .exterMethod + % - The extrapolation method to use with + % scatteredInterpolant. See `help scatteredInterpolant` + % .smoothingLength + % - The smoothing length to use with localSmoothing. + % - AC: can you describe this please + % .fillWith + % - AC: can you describe this please + % .distanceThreshold + % - The distance threshold from known data to truncate the + % interpolated data. + % + % OPENEPDATAINTERPOLATOR Instances of openEpDataInterpolator perform + % spatial interpolation. When instantiated, objects of this type + % describe the interpolation scheme to be used. When ready to perform + % interpolation, the `interpolation(...)` function should be called. + % + % Author: Steven Williams (2021) (Copyright) + % SPDX-License-Identifier: Apache-2.0 + % + % Modifications - + % + % Info on Code Testing: + % --------------------------------------------------------------- + % (1) Perform interpolation of electrogram voltage data using + % localSmoothing with default options. + % + % int = openEpDataInterpolator('localSmoothing'); + % egmX = getElectrogramX(userdata); + % bip = getVoltages(userdata); + % vtx = getVertices(userdata, 'used', false); + % vertexVoltageData = int.interpolate(egmX, bip, vtx); + % + % (2) Visualise the above data + % + % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2]) + % title('Carto voltage data') + % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2], 'data', vertexVoltageData) + % title('OpenEP voltage data') + % --------------------------------------------------------------- + % + % --------------------------------------------------------------- + % code + % --------------------------------------------------------------- + % Properties properties method - end - - properties (Access = private) interMethod exterMethod distanceThreshold @@ -18,14 +76,14 @@ function int = openEpDataInterpolator(varargin) % specify default values + int.method = 'scatteredinterpolant'; int.interMethod = 'linear'; int.exterMethod = 'linear'; int.distanceThreshold = Inf; int.smoothingLength = 10; int.fillWith = NaN; - int.method = 'scatteredinterpolant'; - + % parse input data and warn if conflicts if nargin==1 int.method = varargin{1}; end @@ -71,7 +129,7 @@ end end - % Interpolation + % Perform Interpolation function d1 = interpolate(int, x0,d0,x1) % remove any data with NaNs @@ -99,7 +157,7 @@ case 'localsmoothing' % error checking - % % % AC: is there any appropriate error checking to do + % % % AC: is there any appropriate error checking to do % here? % interpolation From 2fc6fe44c92ba051e4279649ca108a800ad6cb2f Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 9 Apr 2021 16:10:35 +0100 Subject: [PATCH 04/30] Modified default values for scatteredInterpolant.m --- openEpDataInterpolator.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m index 2bc8a00..cd36f63 100644 --- a/openEpDataInterpolator.m +++ b/openEpDataInterpolator.m @@ -50,7 +50,8 @@ % vertexVoltageData = int.interpolate(egmX, bip, vtx); % % (2) Visualise the above data - % + % + % figure;histogram(vertexVoltageData); % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2]) % title('Carto voltage data') % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2], 'data', vertexVoltageData) @@ -78,7 +79,7 @@ % specify default values int.method = 'scatteredinterpolant'; int.interMethod = 'linear'; - int.exterMethod = 'linear'; + int.exterMethod = 'nearest'; int.distanceThreshold = Inf; int.smoothingLength = 10; int.fillWith = NaN; From 15e5f87aadd0f06f3eb9cc0bdd2be2a992c963cd Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 9 Apr 2021 16:44:12 +0100 Subject: [PATCH 05/30] Added blank functions for each of the conduction velocity calculation techniques --- createInterpolator.m | 46 --------------------------------- doCvMapping_CosineFit.m | 48 ++++++++++++++++++++++++++++++++++ doCvMapping_Eikonal.m | 46 +++++++++++++++++++++++++++++++++ doCvMapping_Gradient.m | 51 +++++++++++++++++++++++++++++++++++++ doCvMapping_Omnipole.m | 48 ++++++++++++++++++++++++++++++++++ doCvMapping_RadialBasis.m | 51 +++++++++++++++++++++++++++++++++++++ doCvMapping_Triangulation.m | 48 ++++++++++++++++++++++++++++++++++ getConductionVelocity.m | 14 +++------- 8 files changed, 296 insertions(+), 56 deletions(-) delete mode 100644 createInterpolator.m create mode 100644 doCvMapping_CosineFit.m create mode 100644 doCvMapping_Eikonal.m create mode 100644 doCvMapping_Gradient.m create mode 100644 doCvMapping_Omnipole.m create mode 100644 doCvMapping_RadialBasis.m create mode 100644 doCvMapping_Triangulation.m diff --git a/createInterpolator.m b/createInterpolator.m deleted file mode 100644 index dd99f8b..0000000 --- a/createInterpolator.m +++ /dev/null @@ -1,46 +0,0 @@ -function int = createInterpolator( type, options ) -% CREATEINTERPOLATOR Creates an interpolation scheme for use with OpenEP -% data -% -% Usage: -% int = createInterpolator( type ) -% -% Where: -% type - the type of interpolator required which may be one of, -% `'scatteredInterpolant'`, `'radialbasis'` or `'localsmoothing'` -% int - a function handle to the interpolator -% options - a structure with the following fields: -% options.method -% options.extrapolationmethod -% options.distancethreshold -% -% CREATEINTERPOLATOR creates and interpolation scheme for use with OpenEP -% data -% -% Author: Steven Williams (2021) (Copyright) -% SPDX-License-Identifier: Apache-2.0 -% -% Modifications - -% -% Info on Code Testing: -% --------------------------------------------------------------- -% cvdata = getConductionVelocity( userdata ); -% --------------------------------------------------------------- -% -% --------------------------------------------------------------- -% code -% --------------------------------------------------------------- - -switch(type) - - case 'scatteredInterpolant' - - method = options.method; - extrapolationMethod = options.extrapolationMethod; - distanceThreshold = options.distanceThreshold; - - case 'radialbasis' - - case 'localsmoothing' - -end \ No newline at end of file diff --git a/doCvMapping_CosineFit.m b/doCvMapping_CosineFit.m new file mode 100644 index 0000000..23155fd --- /dev/null +++ b/doCvMapping_CosineFit.m @@ -0,0 +1,48 @@ +function [cv, cvX, interpCv] = doCvMapping_CosineFit( userdata, int ) +% DOCVMAPPING_COSINEFIT Calculates conduction velocities using +% triangulation of electrode activation times +% +% Usage: +% [cv, cvX, interpCv] = doCvMapping_CosineFit( userdata, int ) +% Where: +% userdata - see importcarto_mem +% int - see openEpDataInterpolator.m +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% DOCVMAPPING_COSINEFIT Calculatess conduction velocities using +% cosine fit technique +% +% Author: +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% -----------------------------------s---------------------------- +% +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% first perform global interpolation using cosine fit to +% calculate conduction velocities + +% TODO + +% accept only those conduction velocity values in proximity to electrodes? +% TODO. We should in some way limit cv and cvX only to values that are +% likely to be real; i.e. in close proximity; or at; mapping points. + +% now do interpolation, using the interpolator specified, so we have a full +% dataset at each of the mesh nodes. +vtx = getVertices(userdata, 'used', false); +interpCv = int.interpolate(X, cv, vtx); + +end \ No newline at end of file diff --git a/doCvMapping_Eikonal.m b/doCvMapping_Eikonal.m new file mode 100644 index 0000000..d522182 --- /dev/null +++ b/doCvMapping_Eikonal.m @@ -0,0 +1,46 @@ +function [cv, cvX, interpCv] = doCvMapping_Eikonal( userdata, int ) +% DOCVMAPPING_EIKONAL Calculates conduction velocities using Eikonal +% technique +% +% Usage: +% [cv, cvX, interpCv] = doCvMapping_Eikonal( userdata, int ) +% Where: +% userdata - see importcarto_mem +% int - see openEpDataInterpolator.m +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% DOCVMAPPING_EIKONAL Calculatess conduction velocities using +% omnipoles +% +% Author: +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% -----------------------------------s---------------------------- +% +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% first calculate conduction velocities everywhere +% TODO + +% accept only those conduction velocity values in proximity to electrodes? +% TODO. We should in some way limit cv and cvX only to values that are +% likely to be real; i.e. in close proximity; or at; mapping points. + +% now do interpolation, using the interpolator specified, so we have a full +% dataset at each of the mesh nodes. +vtx = getVertices(userdata, 'used', false); +interpCv = int.interpolate(X, cv, vtx); + +end \ No newline at end of file diff --git a/doCvMapping_Gradient.m b/doCvMapping_Gradient.m new file mode 100644 index 0000000..c26182f --- /dev/null +++ b/doCvMapping_Gradient.m @@ -0,0 +1,51 @@ +function [cv, cvX, interpCv] = doCvMapping_Gradient( userdata, int ) +% DOCVMAPPING_GRADIENT Calculates conduction velocities using +% the gradient of the scalar activation time field +% +% Usage: +% [cv, cvX, interpCv] = doCvMapping_Gradient( userdata, int ) +% Where: +% userdata - see importcarto_mem +% int - see openEpDataInterpolator.m +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% DOCVMAPPING_GRADIENT Calculatess conduction velocities using +% omnipoles +% +% Author: +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% -----------------------------------s---------------------------- +% +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% first perform global interpolation to calculate the activation field +% ?use int (instance of openEpDataInterpolator to do this) +% then we also have several options for how to do this step + +% calculate the gradeint of this field to return the conduction velocities + +% TODO + +% accept only those conduction velocity values in proximity to electrodes? +% TODO. We should in some way limit cv and cvX only to values that are +% likely to be real; i.e. in close proximity; or at; mapping points. + +% now do interpolation, using the interpolator specified, so we have a full +% dataset at each of the mesh nodes. +vtx = getVertices(userdata, 'used', false); +interpCv = int.interpolate(X, cv, vtx); + +end \ No newline at end of file diff --git a/doCvMapping_Omnipole.m b/doCvMapping_Omnipole.m new file mode 100644 index 0000000..93039d7 --- /dev/null +++ b/doCvMapping_Omnipole.m @@ -0,0 +1,48 @@ +function [cv, cvX, interpCv] = doCvMapping_Omnipole( userdata, int ) +% DOCVMAPPING_OMNIPOLE Calculates conduction velocities using +% omnipoles +% +% Usage: +% [cv, cvX, interpCv] = doCvMapping_Omnipole( userdata, int ) +% Where: +% userdata - see importcarto_mem +% int - see openEpDataInterpolator.m +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% DOCVMAPPING_OMNIPOLE Calculatess conduction velocities using +% omnipoles +% +% Author: +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% -----------------------------------s---------------------------- +% +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% first perform global interpolation using omnipoles to +% calculate conduction velocities + +% TODO + +% accept only those conduction velocity values in proximity to electrodes? +% TODO. We should in some way limit cv and cvX only to values that are +% likely to be real; i.e. in close proximity; or at; mapping points. + +% now do interpolation, using the interpolator specified, so we have a full +% dataset at each of the mesh nodes. +vtx = getVertices(userdata, 'used', false); +interpCv = int.interpolate(X, cv, vtx); + +end \ No newline at end of file diff --git a/doCvMapping_RadialBasis.m b/doCvMapping_RadialBasis.m new file mode 100644 index 0000000..0a8a5a0 --- /dev/null +++ b/doCvMapping_RadialBasis.m @@ -0,0 +1,51 @@ +function [cv, cvX, interpCv] = doCvMapping_RadialBasis( userdata, int ) +% DOCVMAPPING_RADIALBASIS Calculates conduction velocities using radial +% basis functions +% +% Usage: +% [cv, cvX, interpCv] = doCvMapping_RadialBasis( userdata ) +% Where: +% userdata - see importcarto_mem +% int - see openEpDataInterpolator.m +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% DOCVMAPPING_RADIALBASIS Calculates conduction velocities using radial +% basis functions +% +% Author: Steven Williams (2021) (Copyright) +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% first perform global interpolation using radial basis function to +% calculate conduction velocities +lats = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); +X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); +cvX = userdata.surface.triRep.X; +[~,~,~,~,cvdata] = RBFConductionVelocity(lats', X', cvX'); +cv = cvdata'; + +% accept only those conduction velocity values in proximity to electrodes? +% TODO. We should in some way limit cv and cvX only to values that are +% likely to be real; i.e. in close proximity; or at; mapping points. + +% now do interpolation, using the interpolator specified, so we have a full +% dataset at each of the mesh nodes. +vtx = getVertices(userdata, 'used', false); +interpCv = int.interpolate(X, cv, vtx); + +end \ No newline at end of file diff --git a/doCvMapping_Triangulation.m b/doCvMapping_Triangulation.m new file mode 100644 index 0000000..5ee2f96 --- /dev/null +++ b/doCvMapping_Triangulation.m @@ -0,0 +1,48 @@ +function [cv, cvX, interpCv] = doCvMapping_Triangulation( userdata, int ) +% DOCVMAPPING_TRIANGULATION Calculates conduction velocities using +% triangulation of electrode activation times +% +% Usage: +% [cv, cvX, interpCv] = doCvMapping_Triangulation( userdata, int ) +% Where: +% userdata - see importcarto_mem +% int - see openEpDataInterpolator.m +% cv - the calculated conduction velocity data, in m/s +% cvX - the Cartesian co-ordinates at which conduction velocity data +% has been calculated. size(cvX) = [length(cv), 3]. +% interpCv - conduction velocity data interpolated across the surface of +% the shell. +% size(interpCv) = [length(userdata.surface.triRep.X), 1]. +% +% DOCVMAPPING_TRIANGULATION Calculatess conduction velocities using +% triangulation +% +% Author: +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% -----------------------------------s---------------------------- +% +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% first perform global interpolation using triangulation to +% calculate conduction velocities + +% TODO + +% accept only those conduction velocity values in proximity to electrodes? +% TODO. We should in some way limit cv and cvX only to values that are +% likely to be real; i.e. in close proximity; or at; mapping points. + +% now do interpolation, using the interpolator specified, so we have a full +% dataset at each of the mesh nodes. +vtx = getVertices(userdata, 'used', false); +interpCv = int.interpolate(X, cv, vtx); + +end \ No newline at end of file diff --git a/getConductionVelocity.m b/getConductionVelocity.m index cfd63ab..d06a896 100644 --- a/getConductionVelocity.m +++ b/getConductionVelocity.m @@ -18,8 +18,8 @@ % `'method'` `'triangulation'` | `'cosinefit'` | `{'radialbasis'}` | % `'omnipole'` | `'gradient'` | `'eikonal'` % `'interpolator'` `{'scatteredInterpolant'}` | `'radialbasis'` | -% `'localsmoothing'` - +% `'localsmoothing'`; or an instance of openEpDataInterpolator +% % GETCONDUCTIONVELOCITY Returns the conduction velocity data of the chamber. % Five methods for calculating conduction velocity are provided as % described below. The conduction velocity data interpolated over the @@ -90,15 +90,10 @@ [cv, cvX, interpCv] = doCvMapping_Triangulation(userdata, int); case 'cosinefit' - [cv, cvX, interpCv] = doCvMapping_CosineFit(userdata); + [cv, cvX, interpCv] = doCvMapping_CosineFit(userdata, int); case 'radialbasis' - [cv, cvX, interpCv] = doCvMapping_RadialBasis(userdata); - - % lats = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); - % X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); - % [~,~,~,~,cvdata] = RBFConductionVelocity(lats', X', userdata.surface.triRep.X'); - % cvdata = cvdata'; + [cv, cvX, interpCv] = doCvMapping_RadialBasis(userdata, int); case 'omnipole' [cv, cvX, interpCv] = doCvMapping_Omnipole(userdata, int); @@ -108,7 +103,6 @@ case 'eikonal' [cv, cvX, interpCv] = doCvMapping_Eikonal(userdata, int); - end end \ No newline at end of file From a563386d0766eb21f4e565253b8763aceac0471a Mon Sep 17 00:00:00 2001 From: Adam Connolly Date: Wed, 14 Apr 2021 08:38:03 +0100 Subject: [PATCH 06/30] updates --- openEpDataInterpolator.m | 42 +++++++++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m index cd36f63..d286e5e 100644 --- a/openEpDataInterpolator.m +++ b/openEpDataInterpolator.m @@ -7,6 +7,10 @@ % int = openEpDataInterpolator(method) % int = openEpDataInterpolator(method, options) % + % Usage (interpolator): + % interpolated_function_at_query_points = int.interpolate(points, + % function_at_points, query_points); + % % Where: % method - the method to use for interpolation, can be % 'scatteredInterpolant', 'localSmoothing' or 'radialBasis' @@ -21,9 +25,14 @@ % scatteredInterpolant. See `help scatteredInterpolant` % .smoothingLength % - The smoothing length to use with localSmoothing. - % - AC: can you describe this please + % Consider a value in the range 5-10. Large values may overly + % smooth the resulting field, and small values may not + % provide enough coverage. % .fillWith - % - AC: can you describe this please + % - The value to assign to the field at query points which + % fall outside the smoothingLength radius. Possible values + % are "nearest" or NaN. "nearest" performs nearest neighbour + % interpolation for query points outside the smoothingLength. % .distanceThreshold % - The distance threshold from known data to truncate the % interpolated data. @@ -70,6 +79,7 @@ distanceThreshold smoothingLength fillWith + rbfConstant end methods @@ -82,7 +92,8 @@ int.exterMethod = 'nearest'; int.distanceThreshold = Inf; int.smoothingLength = 10; - int.fillWith = NaN; + int.fillWith = 'nearest'; + int.rbfConstant = 1; % parse input data and warn if conflicts if nargin==1 @@ -127,6 +138,14 @@ end end end + if isfield(options, 'rbfConstant') + if ~isempty(options.fillWith) + if ~strcmpi(int.method, 'radialBasis') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting rbfConstant property has no effect for interpolation method ' int.method]) + int.rbfConstant = options.rbfConstant; + end + end + end end end @@ -157,18 +176,23 @@ case 'localsmoothing' - % error checking - % % % AC: is there any appropriate error checking to do - % here? - + if isempty(x0) || isempty(x1) || isempty(d0) + error('OPENEP/OPENEPDATAINTERPOLATOR: Missing query or data points!'); + end % interpolation - d1 = localSmoothing(x0, d0, x1, int.smoothingLength, int.fillWith); + [d1, ~] = localSmoothing(x0, d0, x1, int.smoothingLength, int.fillWith); case 'radialbasis' - % error checking + if isempty(x0) || isempty(x1) || isempty(d0) + error('OPENEP/OPENEPDATAINTERPOLATOR: Missing query or data points!'); + end % interpolation + op = rbfcreate(x0', d0', 'RBFFunction', 'multiquadric', 'RBFConstant', int.rbfConstant); + rbfcheck(op); %if this returns anything >10^-9 then consider increasing 'RBFConstant' + [d1, ~] = rbfinterp(x1', op); + d1 = d1'; end % remove interpolated data that is too far from real data From 5f12d161db4c692f6bdfdf2591c11c6e38307d6b Mon Sep 17 00:00:00 2001 From: Iain Sim Date: Fri, 30 Apr 2021 12:05:56 +0100 Subject: [PATCH 07/30] updates to data interpolator class def --- openEpDataInterpolator.m | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m index d286e5e..8fa38a9 100644 --- a/openEpDataInterpolator.m +++ b/openEpDataInterpolator.m @@ -8,7 +8,7 @@ % int = openEpDataInterpolator(method, options) % % Usage (interpolator): - % interpolated_function_at_query_points = int.interpolate(points, + % interpolated_function_at_query_points = int.interpolate(points, % function_at_points, query_points); % % Where: @@ -27,7 +27,7 @@ % - The smoothing length to use with localSmoothing. % Consider a value in the range 5-10. Large values may overly % smooth the resulting field, and small values may not - % provide enough coverage. + % provide enough coverage. % .fillWith % - The value to assign to the field at query points which % fall outside the smoothingLength radius. Possible values @@ -59,7 +59,7 @@ % vertexVoltageData = int.interpolate(egmX, bip, vtx); % % (2) Visualise the above data - % + % % figure;histogram(vertexVoltageData); % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2]) % title('Carto voltage data') @@ -96,25 +96,23 @@ int.rbfConstant = 1; % parse input data and warn if conflicts - if nargin==1 - int.method = varargin{1}; - end + int.method = varargin{1}; if nargin==2 options = varargin{2}; if isfield(options, 'interMethod') if ~isempty(options.interMethod) if ~strcmpi(int.method, 'scatteredInterpolant') warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting interMethod property has no effect for interpolation method ' int.method]); - int.interMethod = options.interMethod; end + int.interMethod = options.interMethod; end end if isfield(options, 'exterMethod') if ~isempty(options.exterMethod) if ~strcmpi(int.method, 'scatteredInterpolant') warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting exterMethod property has no effect for interpolation method ' int.method]); - int.exterMethod = options.exterMethod; end + int.exterMethod = options.exterMethod; end end if isfield(options, 'distanceThreshold') @@ -126,24 +124,24 @@ if ~isempty(options.smoothingLength) if ~strcmpi(int.method, 'localSmoothing') warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) - int.smoothingLength = options.smoothingLength; end + int.smoothingLength = options.smoothingLength; end end if isfield(options, 'fillWith') if ~isempty(options.fillWith) if ~strcmpi(int.method, 'localSmoothing') warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) - int.fillWith = options.fillWith; end + int.fillWith = options.fillWith; end end if isfield(options, 'rbfConstant') if ~isempty(options.fillWith) if ~strcmpi(int.method, 'radialBasis') warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting rbfConstant property has no effect for interpolation method ' int.method]) - int.rbfConstant = options.rbfConstant; end + int.rbfConstant = options.rbfConstant; end end end From ccd23790cd433a8df52ab4f33fedc856e017d7da Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 2 Jul 2021 10:57:14 +0100 Subject: [PATCH 08/30] Added interpolation architecture diagram.drawio --- Untitled Diagram.drawio | 1 + 1 file changed, 1 insertion(+) create mode 100644 Untitled Diagram.drawio diff --git a/Untitled Diagram.drawio b/Untitled Diagram.drawio new file mode 100644 index 0000000..25b5f64 --- /dev/null +++ b/Untitled Diagram.drawio @@ -0,0 +1 @@ +7Zvfb9o6FMf/GqTtoRVJSAiPBdrdSa1UXXa19WkyiSHWTJwaQ6F//T2Ond8JhY4IHjJVW3xiO/E5H5987Xo9a7LafeMoCp6Yj2nP7Pu7njXtmaZhuA78Iy17ZRma2rDkxNeVMsOMvGNt7Gvrhvh4XagoGKOCREWjx8IQe6JgQ5yzt2K1BaPFp0ZoiSuGmYdo1fqT+CJQVtccZvZ/MFkGyZMNZ6TurFBSWY9kHSCfveVM1n3PmnDGhLpa7SaYSuclflHtHhrupi/GcSiOafC+fX33Xgcv9/9NtmK1eF6O/e83upctohs9YBbh8D6aIoG+hwLziFEkGNdDEPvEL+s3sqIohNJ4wUIx03cMKCNKliFce1h2AIYt5oKAS+/0DcEisHoBof4j2rONfP21QN6fpDQOGCfv0C2iuk+4zYWmAwjK15jJlmDug5XjNdR5TnxilExPaFeo+IjWQhs8RimK1mSeDmOF+JKEYyYEW+lK2lkwHLxrjIKRxhYmBWYrLPgequgGCQ16Ohi2Lr9lcJmutgU5sMyBq6HWQC/TrrOYw4UO+wkImBUEeuZYAoxFwKDhHVzPBCfhEi6+RJs5Jd7XnpzG1g7FTsldyoY+gWiFHv4RgJeDeMLJPqZsI517XB9AC/1N8gDGfSg4i2BCf5xskcBfK5BCmEQMD2d/8IRR2c00ZIpaQmnJlIBL8UI0YruOkAe+eIzrTAeZ5V8dLGli0HZB49keEN/HoUSOCZhUii8JU8RgfHEw7TH8QMwn/Vu7Z8OLT6BsZGX4kdW5mLAQxoJIjBoGeN+wBPg4LpvnfxVWTWeSqz+CczRoiU2rhs1SjCmJY6dinKRo41MBXkGoKM4i+kMGfHpjVKJuVaNu1USYojmmz2xNBGGyf67qliJ/qeAa5pGpp63EM2hIPPVfoC9fVQrwYjdtvHjuNyeQLHfgpKVKOydlnbSpyjAfZL2o/KSJHA1Ejqsc0KWncxPsHkdwivrZGbYb9BMoSAEYYT9FGMbZaajzaCh7WBRRA7cmk1m1HLT1oXIaclmcTJ5qlFRzIoGkUN+kyx4nZg/7ZHFTB43TFjPDTty0FtzRkRmhLW3jHsoHWpbslEBY/NYXryX10M338yJRIxbqkEjywtmZGDVohUcGcZmtGBNBXZ7vVMLnVML1iYRkT7OSFdZJ8B9xuJQvl98vaVYKctb/JOX6Xdo4MW2MrlsmGNU92k4nnCu6l9YJRtPuaycULsfEpYWCUd327JRCi0rBMq9PKjRtjvL54mETeiqZHrehwBfzeArJLahOKfzVdqR15VKhuh3ZSYWzhffiWuHgJmOnFS4DxcXFwhHbiDj07+RZGChJIQH+AMsDocmXHEr6Ow56FL7jYpXe2RHxK3f9IkMDPlelafLpjgv7pBDCuH7lC6qVaSflrF1cShqq98Z+5URO6UMOY2Mb7uEDXtE+AIWyxIeiq87xVKObC6fh6HByDDOMbIsvVxdP3d2zxDdTGeUNCdso/VJKDUq3yqiodDQYmLd9Z5T9GRT6NVy4bRiV28ljlE8qjwFA0D5XTU++Mp6pj/7iwEl1zXO/E8Dg+hC4c8qkpMyz6hRA7pfITSnMgfeSv9k+haMjKbQO5xij+uE5F5RDqwilO/wslKNiR1aZ7gbsso6SimyxWONCnXOBZ1VFdbvgXYq7ZAF5xeCZzujW6A+Hg4H+u0BPStOpGNp2CcPyb/6vAcOqSm8Xw+GlOLSvHsOB4972XfgUO5ZrjoauM2yHQ7OcV6+Bw7r1RAnAJWebqNcknvWhay3Ze+nh1hO2Xsp+cmpiXXukp9/WkR6zRlBPzN7dOF5mLRAQD8Wx3Gehcmkx53C1VIsMbUkMdcdVVRV4s3nWrOT03CxudP2H6xYV3cPKts6xdlt+teoOcypvwPowLLjAed3IY/FjuTy90YtOuYqN153p3cyBp66C0yCoJ18iDvuiu1sICxSz/12gkkb2fzSs+/8B \ No newline at end of file From d30c90c73a61014c9a4c8ad1cb65dc40d8c55405 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 2 Jul 2021 13:15:27 +0100 Subject: [PATCH 09/30] created UML folder and moved draw.io diagram --- Untitled Diagram.drawio => UML/interpolation architecture.drawio | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename Untitled Diagram.drawio => UML/interpolation architecture.drawio (100%) diff --git a/Untitled Diagram.drawio b/UML/interpolation architecture.drawio similarity index 100% rename from Untitled Diagram.drawio rename to UML/interpolation architecture.drawio From 75af0e3a992c328aacd5c345d3593c6de505ff2f Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 2 Jul 2021 13:57:29 +0100 Subject: [PATCH 10/30] Update interpolation architecture.drawio --- UML/interpolation architecture.drawio | 115 +++++++++++++++++++++++++- 1 file changed, 114 insertions(+), 1 deletion(-) diff --git a/UML/interpolation architecture.drawio b/UML/interpolation architecture.drawio index 25b5f64..a55074e 100644 --- a/UML/interpolation architecture.drawio +++ b/UML/interpolation architecture.drawio @@ -1 +1,114 @@ -7Zvfb9o6FMf/GqTtoRVJSAiPBdrdSa1UXXa19WkyiSHWTJwaQ6F//T2Ond8JhY4IHjJVW3xiO/E5H5987Xo9a7LafeMoCp6Yj2nP7Pu7njXtmaZhuA78Iy17ZRma2rDkxNeVMsOMvGNt7Gvrhvh4XagoGKOCREWjx8IQe6JgQ5yzt2K1BaPFp0ZoiSuGmYdo1fqT+CJQVtccZvZ/MFkGyZMNZ6TurFBSWY9kHSCfveVM1n3PmnDGhLpa7SaYSuclflHtHhrupi/GcSiOafC+fX33Xgcv9/9NtmK1eF6O/e83upctohs9YBbh8D6aIoG+hwLziFEkGNdDEPvEL+s3sqIohNJ4wUIx03cMKCNKliFce1h2AIYt5oKAS+/0DcEisHoBof4j2rONfP21QN6fpDQOGCfv0C2iuk+4zYWmAwjK15jJlmDug5XjNdR5TnxilExPaFeo+IjWQhs8RimK1mSeDmOF+JKEYyYEW+lK2lkwHLxrjIKRxhYmBWYrLPgequgGCQ16Ohi2Lr9lcJmutgU5sMyBq6HWQC/TrrOYw4UO+wkImBUEeuZYAoxFwKDhHVzPBCfhEi6+RJs5Jd7XnpzG1g7FTsldyoY+gWiFHv4RgJeDeMLJPqZsI517XB9AC/1N8gDGfSg4i2BCf5xskcBfK5BCmEQMD2d/8IRR2c00ZIpaQmnJlIBL8UI0YruOkAe+eIzrTAeZ5V8dLGli0HZB49keEN/HoUSOCZhUii8JU8RgfHEw7TH8QMwn/Vu7Z8OLT6BsZGX4kdW5mLAQxoJIjBoGeN+wBPg4LpvnfxVWTWeSqz+CczRoiU2rhs1SjCmJY6dinKRo41MBXkGoKM4i+kMGfHpjVKJuVaNu1USYojmmz2xNBGGyf67qliJ/qeAa5pGpp63EM2hIPPVfoC9fVQrwYjdtvHjuNyeQLHfgpKVKOydlnbSpyjAfZL2o/KSJHA1Ejqsc0KWncxPsHkdwivrZGbYb9BMoSAEYYT9FGMbZaajzaCh7WBRRA7cmk1m1HLT1oXIaclmcTJ5qlFRzIoGkUN+kyx4nZg/7ZHFTB43TFjPDTty0FtzRkRmhLW3jHsoHWpbslEBY/NYXryX10M338yJRIxbqkEjywtmZGDVohUcGcZmtGBNBXZ7vVMLnVML1iYRkT7OSFdZJ8B9xuJQvl98vaVYKctb/JOX6Xdo4MW2MrlsmGNU92k4nnCu6l9YJRtPuaycULsfEpYWCUd327JRCi0rBMq9PKjRtjvL54mETeiqZHrehwBfzeArJLahOKfzVdqR15VKhuh3ZSYWzhffiWuHgJmOnFS4DxcXFwhHbiDj07+RZGChJIQH+AMsDocmXHEr6Ow56FL7jYpXe2RHxK3f9IkMDPlelafLpjgv7pBDCuH7lC6qVaSflrF1cShqq98Z+5URO6UMOY2Mb7uEDXtE+AIWyxIeiq87xVKObC6fh6HByDDOMbIsvVxdP3d2zxDdTGeUNCdso/VJKDUq3yqiodDQYmLd9Z5T9GRT6NVy4bRiV28ljlE8qjwFA0D5XTU++Mp6pj/7iwEl1zXO/E8Dg+hC4c8qkpMyz6hRA7pfITSnMgfeSv9k+haMjKbQO5xij+uE5F5RDqwilO/wslKNiR1aZ7gbsso6SimyxWONCnXOBZ1VFdbvgXYq7ZAF5xeCZzujW6A+Hg4H+u0BPStOpGNp2CcPyb/6vAcOqSm8Xw+GlOLSvHsOB4972XfgUO5ZrjoauM2yHQ7OcV6+Bw7r1RAnAJWebqNcknvWhay3Ze+nh1hO2Xsp+cmpiXXukp9/WkR6zRlBPzN7dOF5mLRAQD8Wx3Gehcmkx53C1VIsMbUkMdcdVVRV4s3nWrOT03CxudP2H6xYV3cPKts6xdlt+teoOcypvwPowLLjAed3IY/FjuTy90YtOuYqN153p3cyBp66C0yCoJ18iDvuiu1sICxSz/12gkkb2fzSs+/8B \ No newline at end of file + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 63d354a614964c4897196519eac21fbc82a9282a Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 2 Jul 2021 15:53:14 +0100 Subject: [PATCH 11/30] Update interpolation architecture.drawio --- UML/interpolation architecture.drawio | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/UML/interpolation architecture.drawio b/UML/interpolation architecture.drawio index a55074e..bfa6d84 100644 --- a/UML/interpolation architecture.drawio +++ b/UML/interpolation architecture.drawio @@ -1,6 +1,6 @@ - + - + From 0bdf1136400863645805c5b269b28537cc0db7d3 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Sun, 18 Jul 2021 18:36:40 +0100 Subject: [PATCH 12/30] Finished first beta of conduction velocity architecture --- doCvMapping_Gradient.m | 50 +++++--- doCvMapping_RadialBasis.m | 5 +- getConductionVelocity.m | 2 +- getVerticesNearMappingPoints.m | 67 +++++++++++ openEpDataInterpolator.m | 154 +++++++----------------- openEpDataInterpolator_bu.m | 206 +++++++++++++++++++++++++++++++++ openEpScatteredInterpolant.m | 53 +++++++++ private/trigrad.m | 100 ++++++++++++++++ 8 files changed, 508 insertions(+), 129 deletions(-) create mode 100644 getVerticesNearMappingPoints.m create mode 100644 openEpDataInterpolator_bu.m create mode 100644 openEpScatteredInterpolant.m create mode 100644 private/trigrad.m diff --git a/doCvMapping_Gradient.m b/doCvMapping_Gradient.m index c26182f..1b9ed11 100644 --- a/doCvMapping_Gradient.m +++ b/doCvMapping_Gradient.m @@ -24,28 +24,52 @@ % % Info on Code Testing: % -----------------------------------s---------------------------- -% +% load openep_dataset_1.mat +% int = openEpDataInterpolator +% [cv, cvX, interpCv] = getConductionVelocity(userdata, 'method', 'gradient', 'interpolator', int); +% histogram(cv) +% drawMap(userdata, 'data', interpCv, 'type', 'cv') % --------------------------------------------------------------- % % --------------------------------------------------------------- % code % --------------------------------------------------------------- -% first perform global interpolation to calculate the activation field -% ?use int (instance of openEpDataInterpolator to do this) -% then we also have several options for how to do this step +CVLIMIT = 10; %m/s +DISTANCETHRESHOLD = 10; %mm + +% get relevant mapping points +X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); + +% calculate the local activation times at each mapping point +mapAnnot = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); +refAnnot = userdata.electric.annotations.referenceAnnot(getMappingPointsWithinWoI(userdata)); +lats = mapAnnot - refAnnot; + +% perform global interpolation of the activation field +latMap = int.interpolate(X, lats, getVertices(userdata)); + +% calculate the gradeint of this field +[cvX, grad] = trigrad(getMesh(userdata), latMap); -% calculate the gradeint of this field to return the conduction velocities +% calculate conduction velocities +sgrad = grad .* grad; +dp = sum(sgrad,2); +cv = 1./sqrt(dp); -% TODO +% accept only those conduction velocity values in proximity to electrodes +vtx = getVerticesNearMappingPoints(userdata, DISTANCETHRESHOLD); +cv(~vtx) = []; +cvX(~vtx,:) = []; +disp(['OPENEP/DOCVMAPPING_GRADIENT: ' num2str(sum(~vtx)) ' CV values were removed which were more than ' num2str(DISTANCETHRESHOLD) 'mm from a mapping point']); -% accept only those conduction velocity values in proximity to electrodes? -% TODO. We should in some way limit cv and cvX only to values that are -% likely to be real; i.e. in close proximity; or at; mapping points. +% then remove any non physiological values over the CVLIMIT +isOverCvLimit = cv>CVLIMIT; +cv(isOverCvLimit) = []; +cvX(isOverCvLimit,:) = []; +disp(['OPENEP/DOCVMAPPING_GRADIENT: ' num2str(sum(isOverCvLimit)) ' CV values were removed which were greater than ' num2str(CVLIMIT) 'm/s']); -% now do interpolation, using the interpolator specified, so we have a full -% dataset at each of the mesh nodes. -vtx = getVertices(userdata, 'used', false); -interpCv = int.interpolate(X, cv, vtx); +% interpolate +interpCv = int.interpolate(cvX, cv, getVertices(userdata)); end \ No newline at end of file diff --git a/doCvMapping_RadialBasis.m b/doCvMapping_RadialBasis.m index 0a8a5a0..1a8079e 100644 --- a/doCvMapping_RadialBasis.m +++ b/doCvMapping_RadialBasis.m @@ -33,8 +33,11 @@ % first perform global interpolation using radial basis function to % calculate conduction velocities -lats = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); +mapAnnot = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); +refAnnot = userdata.electric.annotations.referenceAnnot(getMappingPointsWithinWoI(userdata)); +lats = mapAnnot - refAnnot; X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); + cvX = userdata.surface.triRep.X; [~,~,~,~,cvdata] = RBFConductionVelocity(lats', X', cvX'); cv = cvdata'; diff --git a/getConductionVelocity.m b/getConductionVelocity.m index d06a896..ca1e419 100644 --- a/getConductionVelocity.m +++ b/getConductionVelocity.m @@ -80,7 +80,7 @@ % first create an interpolator if ischar(interpolator) - int = createInterpolator(interpolator); + int = openEpDataInterpolator(interpolator); else int = interpolator; end diff --git a/getVerticesNearMappingPoints.m b/getVerticesNearMappingPoints.m new file mode 100644 index 0000000..7f29b96 --- /dev/null +++ b/getVerticesNearMappingPoints.m @@ -0,0 +1,67 @@ +function [ vtx ] = getVerticesNearMappingPoints( userdata, distanceThreshold, varargin ) +% GETVERTICESNEARMAPPINGPOINTS identifies mesh vertices which are within +% distanceThreshold of an actual mapping point +% +% Usage: +% [ vol ] = getVerticesNearMappingPoints( userdata, distanceThreshold ) +% Where: +% userdata - see importcarto_mem +% distanceThreshold - the distance threshold to be applied +% vtx - logical array of size length(userdata.surface.triRep.X) which +% identifies the vertices within distanceThreshold of an actual +% mapping point. +% +% GETVERTICESNEARMAPPINGPOINTS accepts the following parameter-value pairs +% 'method' {'linear'}|'geodesic' TODO: N.B. geodesic not yet implemented +% +% GETVERTICESNEARMAPPINGPOINTS Can be used to identify mesh vertices which +% are within a particular distanceThreshold of an original mapping point. +% For example, the function can identify all the mesh vertices within 10mm +% of a local activation time mapping point. This is useful in order to +% filter interpolated data by only those points that are close to points of +% actual data. Two methods are provided for calculating distances: +% Euclidean distances and geodesic distances. Note that geodesic distances +% is not yet implemented. +% +% Author: Steven Williams (2021) +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% test code +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +nStandardArgs = 2; % UPDATE VALUE +method = 'linear'; +if nargin > nStandardArgs + for i = 1:2:nargin-nStandardArgs + switch varargin{i} + case 'method' + method = varargin{i+1}; + end + end +end + +% identify the vertices +% x1 are interpolated data i.e. all the vertices +% x0 real data i.e. co-ordinates of the mapping points +x1 = getVertices(userdata); +x0 = userdata.electric.egmSurfX; + +if strcmpi(method, 'geodesic') + error('OPENEP/GETVERTICESNEARMAPPINGPOINTS: Geodesic distances not yet implemented'); + %TODO: Geodesic functionality needs to be added to distBetweenPoints.m; when working, this function will work too +end + +id = knnsearch(x0, x1); +cPts = x0(id,:); %c for closest +d = distBetweenPoints(cPts, x1, 'method', method); +vtx = ones(size(d)); +vtx(d>distanceThreshold) = 0; +vtx = logical(vtx); + +end \ No newline at end of file diff --git a/openEpDataInterpolator.m b/openEpDataInterpolator.m index 8fa38a9..7c09565 100644 --- a/openEpDataInterpolator.m +++ b/openEpDataInterpolator.m @@ -70,137 +70,63 @@ % --------------------------------------------------------------- % code % --------------------------------------------------------------- - + % Properties properties method - interMethod - exterMethod distanceThreshold - smoothingLength - fillWith - rbfConstant + end + + properties (Access = private) + call_interpolator end methods % Contructor - function int = openEpDataInterpolator(varargin) + function obj = openEpDataInterpolator(method) + + possibleMethods = {'scatteredinterpolant' ... + , 'localsmoothing' ... + , 'radialbasis' ... + }; - % specify default values - int.method = 'scatteredinterpolant'; - int.interMethod = 'linear'; - int.exterMethod = 'nearest'; - int.distanceThreshold = Inf; - int.smoothingLength = 10; - int.fillWith = 'nearest'; - int.rbfConstant = 1; + defaultMethod = 'scatteredinterpolant'; - % parse input data and warn if conflicts - int.method = varargin{1}; - if nargin==2 - options = varargin{2}; - if isfield(options, 'interMethod') - if ~isempty(options.interMethod) - if ~strcmpi(int.method, 'scatteredInterpolant') - warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting interMethod property has no effect for interpolation method ' int.method]); - end - int.interMethod = options.interMethod; - end - end - if isfield(options, 'exterMethod') - if ~isempty(options.exterMethod) - if ~strcmpi(int.method, 'scatteredInterpolant') - warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting exterMethod property has no effect for interpolation method ' int.method]); - end - int.exterMethod = options.exterMethod; - end - end - if isfield(options, 'distanceThreshold') - if ~isempty(options.distanceThreshold) - int.distanceThreshold = options.distanceThreshold; - end - end - if isfield(options, 'smoothingLength') - if ~isempty(options.smoothingLength) - if ~strcmpi(int.method, 'localSmoothing') - warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) - end - int.smoothingLength = options.smoothingLength; - end - end - if isfield(options, 'fillWith') - if ~isempty(options.fillWith) - if ~strcmpi(int.method, 'localSmoothing') - warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) - end - int.fillWith = options.fillWith; - end - end - if isfield(options, 'rbfConstant') - if ~isempty(options.fillWith) - if ~strcmpi(int.method, 'radialBasis') - warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting rbfConstant property has no effect for interpolation method ' int.method]) + obj.distanceThreshold = 10; + + if nargin == 0 + obj.method = defaultMethod; + obj.call_interpolator = openEpScatteredInterpolant; + else + if nargin > 1 + error('Too many input arguments!'); + else + if ~ismember(lower(method), possibleMethods) + error('Choose from ''scatteredinterpolant'', ''localsmoothing'' or ''radialbasis'''); + else + obj.method = lower(method); + switch obj.method + case 'scatteredinterpolant' + obj.call_interpolator = openEpScatteredInterpolant; + case 'localsmoothing' + obj.call_interpolator = openEpLocalSmoothingInterpolant; + case 'radialbasis' + obj.call_interpolator = openEpRadialBasisInterpolant; end - int.rbfConstant = options.rbfConstant; end - end + end end + end - % Perform Interpolation - function d1 = interpolate(int, x0,d0,x1) + function f_q = interpolate(obj, x, f_x, q) - % remove any data with NaNs - tempData = d0; - tempCoords = x0; - iNaN = isnan(tempData); - tempData(iNaN) = []; - tempCoords(iNaN,:) = []; - d0 = tempData; - x0 = tempCoords; - - % perform interpolation - switch lower(int.method) - case 'scatteredinterpolant' - - % error checking - if numel(d0) < 4 - error('OPENEP/OPENEPDATAINTERPOLATOR: Cannot perform interpolation with less than four data points'); - end - - % interpolation - F = scatteredInterpolant(x0(:,1), x0(:,2), x0(:,3), d0, int.interMethod, int.exterMethod); - d1 = F(x1); - - case 'localsmoothing' - - if isempty(x0) || isempty(x1) || isempty(d0) - error('OPENEP/OPENEPDATAINTERPOLATOR: Missing query or data points!'); - end - % interpolation - [d1, ~] = localSmoothing(x0, d0, x1, int.smoothingLength, int.fillWith); - - case 'radialbasis' - - if isempty(x0) || isempty(x1) || isempty(d0) - error('OPENEP/OPENEPDATAINTERPOLATOR: Missing query or data points!'); - end - - % interpolation - op = rbfcreate(x0', d0', 'RBFFunction', 'multiquadric', 'RBFConstant', int.rbfConstant); - rbfcheck(op); %if this returns anything >10^-9 then consider increasing 'RBFConstant' - [d1, ~] = rbfinterp(x1', op); - d1 = d1'; - end + f_q = obj.call_interpolator.interpolate(x, f_x, q); - % remove interpolated data that is too far from real data - id = knnsearch(x0, x1); - cPts = x0(id,:); %c for closest - d = distBetweenPoints(cPts, x1); - thresholdDistance = zeros(size(d)); - thresholdDistance(d>int.distanceThreshold) = 1; - nanset = logical(thresholdDistance); - d1(nanset) = NaN; + end + + function disp(obj) + display(obj.call_interpolator); end end end \ No newline at end of file diff --git a/openEpDataInterpolator_bu.m b/openEpDataInterpolator_bu.m new file mode 100644 index 0000000..2b53048 --- /dev/null +++ b/openEpDataInterpolator_bu.m @@ -0,0 +1,206 @@ +classdef openEpDataInterpolator_bu + % OPENEPDATAINTERPOLATOR Creates objects for performing spatial + % interpolation for OpenEP data + % + % Usage (Constructor): + % int = openEpDataInterpolator() + % int = openEpDataInterpolator(method) + % int = openEpDataInterpolator(method, options) + % + % Usage (interpolator): + % interpolated_function_at_query_points = int.interpolate(points, + % function_at_points, query_points); + % + % Where: + % method - the method to use for interpolation, can be + % 'scatteredInterpolant', 'localSmoothing' or 'radialBasis' + % options - a structure containing options pertaining to each + % interpolation method. Options which are not relevant to a + % particular interpolation method will be ignored. + % .interMethod + % - The interpolation method to use with + % scatteredInterpolant. See `help scatteredInterpolant` + % .exterMethod + % - The extrapolation method to use with + % scatteredInterpolant. See `help scatteredInterpolant` + % .smoothingLength + % - The smoothing length to use with localSmoothing. + % Consider a value in the range 5-10. Large values may overly + % smooth the resulting field, and small values may not + % provide enough coverage. + % .fillWith + % - The value to assign to the field at query points which + % fall outside the smoothingLength radius. Possible values + % are "nearest" or NaN. "nearest" performs nearest neighbour + % interpolation for query points outside the smoothingLength. + % .distanceThreshold + % - The distance threshold from known data to truncate the + % interpolated data. + % + % OPENEPDATAINTERPOLATOR Instances of openEpDataInterpolator perform + % spatial interpolation. When instantiated, objects of this type + % describe the interpolation scheme to be used. When ready to perform + % interpolation, the `interpolation(...)` function should be called. + % + % Author: Steven Williams (2021) (Copyright) + % SPDX-License-Identifier: Apache-2.0 + % + % Modifications - + % + % Info on Code Testing: + % --------------------------------------------------------------- + % (1) Perform interpolation of electrogram voltage data using + % localSmoothing with default options. + % + % int = openEpDataInterpolator('localSmoothing'); + % egmX = getElectrogramX(userdata); + % bip = getVoltages(userdata); + % vtx = getVertices(userdata, 'used', false); + % vertexVoltageData = int.interpolate(egmX, bip, vtx); + % + % (2) Visualise the above data + % + % figure;histogram(vertexVoltageData); + % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2]) + % title('Carto voltage data') + % figure;drawMap(userdata, 'type', 'bip', 'coloraxis', [0.05 2], 'data', vertexVoltageData) + % title('OpenEP voltage data') + % --------------------------------------------------------------- + % + % --------------------------------------------------------------- + % code + % --------------------------------------------------------------- + + % Properties + properties + method + interMethod + exterMethod + distanceThreshold + smoothingLength + fillWith + rbfConstant + end + + methods + % Contructor + function int = openEpDataInterpolator(varargin) + + % specify default values + int.method = 'scatteredinterpolant'; + int.interMethod = 'linear'; + int.exterMethod = 'nearest'; + int.distanceThreshold = Inf; + int.smoothingLength = 10; + int.fillWith = 'nearest'; + int.rbfConstant = 1; + + % parse input data and warn if conflicts + int.method = varargin{1}; + if nargin==2 + options = varargin{2}; + if isfield(options, 'interMethod') + if ~isempty(options.interMethod) + if ~strcmpi(int.method, 'scatteredInterpolant') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting interMethod property has no effect for interpolation method ' int.method]); + end + int.interMethod = options.interMethod; + end + end + if isfield(options, 'exterMethod') + if ~isempty(options.exterMethod) + if ~strcmpi(int.method, 'scatteredInterpolant') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting exterMethod property has no effect for interpolation method ' int.method]); + end + int.exterMethod = options.exterMethod; + end + end + if isfield(options, 'distanceThreshold') + if ~isempty(options.distanceThreshold) + int.distanceThreshold = options.distanceThreshold; + end + end + if isfield(options, 'smoothingLength') + if ~isempty(options.smoothingLength) + if ~strcmpi(int.method, 'localSmoothing') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) + end + int.smoothingLength = options.smoothingLength; + end + end + if isfield(options, 'fillWith') + if ~isempty(options.fillWith) + if ~strcmpi(int.method, 'localSmoothing') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting fillWith property has no effect for interpolation method ' int.method]) + end + int.fillWith = options.fillWith; + end + end + if isfield(options, 'rbfConstant') + if ~isempty(options.fillWith) + if ~strcmpi(int.method, 'radialBasis') + warning(['OPENEP/OPENEPDATAINTERPOLATOR: Setting rbfConstant property has no effect for interpolation method ' int.method]) + end + int.rbfConstant = options.rbfConstant; + end + end + end + end + + % Perform Interpolation + function d1 = interpolate(int, x0,d0,x1) + + % remove any data with NaNs + tempData = d0; + tempCoords = x0; + iNaN = isnan(tempData); + tempData(iNaN) = []; + tempCoords(iNaN,:) = []; + d0 = tempData; + x0 = tempCoords; + + % perform interpolation + switch lower(int.method) + case 'scatteredinterpolant' + + % error checking + if numel(d0) < 4 + error('OPENEP/OPENEPDATAINTERPOLATOR: Cannot perform interpolation with less than four data points'); + end + + % interpolation + F = scatteredInterpolant(x0(:,1), x0(:,2), x0(:,3), d0, int.interMethod, int.exterMethod); + d1 = F(x1); + + case 'localsmoothing' + + if isempty(x0) || isempty(x1) || isempty(d0) + error('OPENEP/OPENEPDATAINTERPOLATOR: Missing query or data points!'); + end + % interpolation + [d1, ~] = localSmoothing(x0, d0, x1, int.smoothingLength, int.fillWith); + + case 'radialbasis' + + if isempty(x0) || isempty(x1) || isempty(d0) + error('OPENEP/OPENEPDATAINTERPOLATOR: Missing query or data points!'); + end + + % interpolation + op = rbfcreate(x0', d0', 'RBFFunction', 'multiquadric', 'RBFConstant', int.rbfConstant); + rbfcheck(op); %if this returns anything >10^-9 then consider increasing 'RBFConstant' + [d1, ~] = rbfinterp(x1', op); + d1 = d1'; + end + + % remove interpolated data that is too far from real data + id = knnsearch(x0, x1); + cPts = x0(id,:); %c for closest + d = distBetweenPoints(cPts, x1); + thresholdDistance = zeros(size(d)); + thresholdDistance(d>int.distanceThreshold) = 1; + nanset = logical(thresholdDistance); + d1(nanset) = NaN; + end + end +end \ No newline at end of file diff --git a/openEpScatteredInterpolant.m b/openEpScatteredInterpolant.m new file mode 100644 index 0000000..8d4e9c6 --- /dev/null +++ b/openEpScatteredInterpolant.m @@ -0,0 +1,53 @@ +classdef openEpScatteredInterpolant + + properties + interMethod = 'linear'; + exterMethod = 'nearest'; + distanceThreshold; + end + + methods + % Default implicit constructor. + function obj = openEpScatteredInterpolantd(distanceThreshold) + obj.distanceThreshold = distanceThreshold; + end + + % Main method + function f_q = interpolate(obj, x, f_x, q) + + F = scatteredInterpolant(x(:,1), x(:,2), x(:,3), f_x, obj.interMethod, obj.exterMethod); + f_q = F(q); + + end + end + + % Getter and setter methods + methods + + function obj = set.interMethod(obj, value) + if ismember(value, {'linear', 'nearest', 'natural'}) + obj.interMethod = value; + else + error('Interpolation method must be one of ''linear'', ''nearest'' or ''natural'''); + end + end + + function value = get.interMethod(obj) + value = obj.interMethod; + end + + function obj = set.exterMethod(obj, value) + if ismember(value, {'linear', 'nearest', 'none'}) + obj.exterMethod = value; + else + error('Extrapolation method must be one of ''linear'', ''nearest'' or or ''none'''); + end + end + + function value = get.exterMethod(obj) + value = obj.exterMethod; + end + + end + +end \ No newline at end of file diff --git a/private/trigrad.m b/private/trigrad.m new file mode 100644 index 0000000..1c37243 --- /dev/null +++ b/private/trigrad.m @@ -0,0 +1,100 @@ +function [locations, gradients] = trigrad(tri, sc) +% TRIGRAD calculates the gradient of a scalar field across a TriRep object +% Usage: +% [locations gradients] = trigrad(tri, sc) +% Where: +% locations contains the position at the centroid of each face +% gradients contains the vector of grad(sc), corresponding to locations +% tri is the trirep object containing the geometry +% sc is the value of the scalar field at each point in tri.X +% g is the gradient +% TRIGRAD calculates the gradient of sc across the surface defined by tri. +% The geometry can be a surface in 2d or 3d. + +% Author: Nick Linton (2009) +% Modifications - + +% Info on Code Testing: + % --------------------- + % test code + % --------------------- + % %intiate random things + % if true + % nPoints = 2000; + % p = rand(nPoints,3)*2*pi; + % for i = 1:nPoints + % p(i,:) = [1 1 1] * rotationmatrix(p(i,1), p(i,2), p(i,3)); + % end + % end + % p(:,1) = p(:,1)-1; + % + % tri = convhulln(p); + % tri = TriRep(tri,p); + % + % SENSITIVITY = 0.1; + % + % for i = 1:length(tri.X(:,1)) + % potential = tri.X(:,1) .* exp( -tri.X(:,1).^2 - tri.X(:,2).^2 ); + % end + % [locations gradients] = trigrad(tri, potential); + % + % close all + % figure + % hold on + % colormap(jet(64)) + % c = (potential-min(potential)) ./ (max(potential)-min(potential)) * 63 + 1; + % h = trisurf(tri); + % set(h, 'CData',c ,'CDataMapping','direct', 'FaceColor', 'interp') + % for i = 1:size(locations,1) + % start = locations(i,:); + % loc = [start; start]; + % loc(2,:) = loc(2,:) + SENSITIVITY * gradients(i,:); + % + % plot3(start(:,1), start(:,2), start(:,3), 'Marker','*') + % plot3(loc(:,1), loc(:,2), loc(:,3)) + % end + % axis equal + % axis vis3d + + + [~, locations] = tricentroid(tri); + gradients = zeros(size(tri.Triangulation,1), 3); + x = tri.X; + nDim = size(x,2); + if nDim == 2 + %then we only have 2 dimensions - so expand to 3, and retract later + x = [x zeros(size(x,1),1)]; + end + for i = 1:size(tri.Triangulation,1) + gradients(i,:) = grad(x(tri.Triangulation(i,:),:), sc(tri.Triangulation(i,:))); + end + if nDim == 2 + gradients(:,3) = []; + end +end + +function g = grad(p, s) + % first find two vectors that are orthogonal and on the surface of p1/2/3 + u = p(2,:)-p(1,:); + u = u / norm(u); + %u = reshape(u,1,length(u)); %make sure it's a row vector [it always is] + + v = p(3,:)-p(1,:); + v = v - u.*dot(u,v); %v is the component of v normal to u (u has already been normalised) + v = v / norm(v); + %v = reshape(v,1,length(v)); %make sure it's a row vector [it always is] + + % now express the change in s across the triangle + % ds = a1 * dp.u + a2 * dp.v + a = [ dot(p(2,:)-p(1,:),u) , dot(p(2,:)-p(1,:),v) ; dot(p(3,:)-p(1,:),u) , dot(p(3,:)-p(1,:),v) ] \ [ (s(2)-s(1)) ; (s(3)-s(1)) ]; + % so now, change in activation time is given by + % s-s(1) = [u v]*a + % so the vector [u v]*a is the gradient of s + g = a(1)*u + a(2)*v; +end + + + + + + \ No newline at end of file From 3d8aa63add806f65c3a1ea3d5ea041702aa5da3e Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Mon, 19 Jul 2021 17:29:07 +0100 Subject: [PATCH 13/30] Added triangulation method although this needs converted into "OpenEP" format still --- computeCVtriangulation.m | 45 ++++++++++++++++++++++++++++++++++++++++ doCvMapping_Gradient.m | 2 +- 2 files changed, 46 insertions(+), 1 deletion(-) create mode 100644 computeCVtriangulation.m diff --git a/computeCVtriangulation.m b/computeCVtriangulation.m new file mode 100644 index 0000000..8b53791 --- /dev/null +++ b/computeCVtriangulation.m @@ -0,0 +1,45 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Computes CV magnitude and direction using triangulation method +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [v v_vec] = computeCVtriangulation(p,tp,q,tq,r,tr) +% inputs: +% coordinate vectors of 3 points defining a triangle: p, q, r i.e. p=[px,py,pz] +% local activation timings of each point: tp, tq, tr +% outputs CV as magntidue as well as normalised velocity vector + +% Defines lengths of triange edges +x_pq = [q(1)-p(1);q(2)-p(2);q(3)-p(3)]; +x_pr = [r(1)-p(1);r(2)-p(2);r(3)-p(3)]; +x_qr = [r(1)-q(1);r(2)-q(2);r(3)-q(3)]; + +% Defines relative activation timings +t_pq = abs(tp - tq); +t_pr = abs(tp - tr); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Computes velocity magnitude +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Computes angle theta (between x_pq and x_pr) using cosine rule +theta = (norm(x_pq)^2 + norm(x_pr)^2 - norm(x_qr)^2)/(2*norm(x_pq)*norm(x_pr)); + +% Computes angle alpha (between x_pq and velocity vector) +alpha = arctan( ( (t_pr*norm(x_pq))/(t_pq*norm(x_pr)) - cos(theta) ) / sin(theta) ); + +% Computes velocity +v = norm(x_pq)*cos(theta)/t_pq; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Computes velocity as vector +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Defines normal to plane (as unit vector) +n_plane_full = cross(x_pr,x_pq); +n_plane = n_plane_full/norm(n_plane_full); + +% Defines vector xs (defined as a vector perpendicular to x_pq, within the plane, from p which intersects with the velocity vector) +x_ps = cross(n_plane,x_pq)*tan(alpha); + +% Computes the velocity vector +vec_full = x_pq - x_ps; +vec = vec_full/norm(vec_full); + + diff --git a/doCvMapping_Gradient.m b/doCvMapping_Gradient.m index 1b9ed11..21d3f14 100644 --- a/doCvMapping_Gradient.m +++ b/doCvMapping_Gradient.m @@ -63,7 +63,7 @@ cvX(~vtx,:) = []; disp(['OPENEP/DOCVMAPPING_GRADIENT: ' num2str(sum(~vtx)) ' CV values were removed which were more than ' num2str(DISTANCETHRESHOLD) 'mm from a mapping point']); -% then remove any non physiological values over the CVLIMIT +% remove any non physiological values over the CVLIMIT isOverCvLimit = cv>CVLIMIT; cv(isOverCvLimit) = []; cvX(isOverCvLimit,:) = []; From 05952149530dfe2782b0fd29a6a0d559295b2b61 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Fri, 10 Sep 2021 16:04:54 +0100 Subject: [PATCH 14/30] improved handling of trirep/triangulation. colouring shell based on geodesic distances. plotting tags. --- drawMap.m | 28 +++++++++-- getData2MeshGeodesicDistances.m | 83 +++++++++++++++++++++++++++++++++ getFaces.m | 8 +++- getMesh.m | 22 +++++++-- getNormals.m | 48 +++++++++++++++++++ getVertices.m | 8 +++- plotTag.m | 50 +++++++++++++++----- private/colorShell.m | 13 ++++++ private/drawFreeBoundary.m | 14 ++++-- 9 files changed, 250 insertions(+), 24 deletions(-) create mode 100644 getData2MeshGeodesicDistances.m create mode 100644 getNormals.m diff --git a/drawMap.m b/drawMap.m index 9752ddf..5c88d42 100644 --- a/drawMap.m +++ b/drawMap.m @@ -11,7 +11,7 @@ % DRAWMAP accepts the following parameter-value pairs % 'data' {[]} | [d] % - Where d is a vector of data values and size(d) equals numel(userdata.surface.triRep.X) -% 'type' {'act'} | 'bip' | 'force' | 'uni' | 'none' | 'cv' +% 'type' {'act'} | 'bip' | 'force' | 'uni' | 'none' | 'cv' | 'geodesic' | 'wallthickness' % - Specifies type of map - activation, bipolar or unipolar voltage % 'coloraxis' {[]} | [a b] % - Where a and b are real numbers. See help colorShell @@ -74,8 +74,8 @@ orientation = inputs.orientation; DISTANCETHRESH = inputs.colorfillthreshold; -if ~any(strcmpi(type, {'act', 'bip', 'force', 'uni', 'none', 'cv'})) - error('DRAWMAP: Invalid parameter-value pair; type must be one of act, bip, force, uni, none') +if ~any(strcmpi(type, {'act', 'bip', 'force', 'uni', 'none', 'cv', 'geodesic', 'wallthickness'})) + error('DRAWMAP: Invalid parameter-value pair; type must be one of act, bip, force, uni, none, cv or geodesic') end @@ -124,6 +124,26 @@ else DATA = getConductionVelocity(userdata); end + case 'geodesic' + DATATYPE = 'geodesic'; + if ~isempty(DATAMANUAL) + DATA = DATAMANUAL; + else + % by default draw the geodesics from the first electrode positon + distances = getData2MeshGeodesicDistances(userdata); + DATA = NaN(size(getVertices(userdata),1),10); + int = openEpDataInterpolator(); + for i = 1:20%size(distances,1) + DATA(:,i) = int.interpolate(getVertices(userdata,'used',true), distances(i,:)', getVertices(userdata,'used',false)); + end + end + case 'wallthickness' + DATATYPE = 'wallthickness'; + if ~isempty(DATAMANUAL) + DATA = DATAMANUAL; + else + error('OPENEP/drawMap: automatic calculation of wall thickness data not supported') + end end if ~strcmpi(type, 'none') @@ -148,7 +168,7 @@ ); else disp('no data to color the shell with') - end + end end % Adjust the light/material diff --git a/getData2MeshGeodesicDistances.m b/getData2MeshGeodesicDistances.m new file mode 100644 index 0000000..c0627ae --- /dev/null +++ b/getData2MeshGeodesicDistances.m @@ -0,0 +1,83 @@ +function [distances, isVertUsed, userdata] = getData2MeshGeodesicDistances(userdata, varargin) +% GETDATA2MESHGEODESICDISTANCES Calculates geodesic distances from data +% points to mesh points +% +% Usage: +% distances = getData2meshGeodesicDistances(userdata) +% [distances, isVertUsed] = getData2meshGeodesicDistances(userdata) +% [distances, isVertUsed, userdata] = getData2meshGeodesicDistances(userdata) +% +% Where: +% userdata - an openEP data structure +% distances - an n x m array of distances, where n is the number of +% electrode data points and m is the number of mesh nodes +% isVertUsed - the vertices that are referenced by the triangulation +% +% GETDATA2MESHGEODESICDISTANCES repacks the original mesh to ensure that +% there are no vertices unreferenced by the triangulation. The vertices for +% which geodesic distances have been calculated are returned in isVertUsed. +% If three output argumetns are specified then a modified userdata +% structure containing the computed geodesic distances is also returned. +% +% See also: repack.m +% +% Author: Steven Williams (2016) +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% int = openEpDataInterpolator +% distances = getData2MeshGeodesicDistances(userdata); +% d = int.interpolate(getVertices(userdata,'used',true), distances(1,:)', getVertices(userdata,'used',false)); +% drawMap(userdata, 'type', 'bip', 'coloraxis', [0 120], 'data', d) +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +global geodesic_library + +if isfield(userdata.surface, 'geodesicDistances') + distances = userdata.surface.geodesicDistances.distances; + isVertUsed = userdata.surface.geodesicDistances.isVertUsed; +else + % compute geodesic distances + + tic + surfaceTriangulation = getMesh(userdata, 'type', 'triangulation'); + [surfaceTriangulation, isVertUsed] = repack(surfaceTriangulation); + geodesic_library = 'geodesic_matlab_api'; + mesh = geodesic_new_mesh(surfaceTriangulation.Points, surfaceTriangulation.ConnectivityList); + algorithm = geodesic_new_algorithm(mesh, 'exact'); + + numPts = getNumPts(userdata); + [~, startPts] = getElectrogramX(userdata, 'type', 'bip'); + + distances = NaN(numPts, size(surfaceTriangulation.Points,1)); + for i = 1:numPts + source_points{1} = geodesic_create_surface_point('vertex', i, startPts(i,:)); + geodesic_propagate(algorithm, source_points); + [~, distances(i,:)] = geodesic_distance_and_source(algorithm); + end + + geodesic_delete() + + T = toc; + disp([ num2str(size(distances,2)) ... + , ' geodesic distances from ' ... + , num2str(size(distances,1)) ... + , ' datapoints calculated in ' ... + , num2str(T) ... + , ' seconds.' ... + ]); + + % store the data + userdata.surface.geodesicDistances.distances = distances; + userdata.surface.geodesicDistances.isVertUsed = isVertUsed; + userdata.surface.geodesicDistances.computeTime = T; +end + +clear geodesic_library + +end \ No newline at end of file diff --git a/getFaces.m b/getFaces.m index d9805d8..65775a5 100644 --- a/getFaces.m +++ b/getFaces.m @@ -23,6 +23,12 @@ % code % --------------------------------------------------------------- -faces = userdata.surface.triRep.Triangulation; +if isa(userdata.surface.triRep, 'TriRep') + faces = userdata.surface.triRep.Triangulation; +elseif isa(userdata.surface.triRep, 'triangulation') + faces = userdata.surface.triRep.ConnectivityList; +else + error('OPENEP/getFaces: userdata.surface.TriRep should be a TriRep or a triangulation object') +end end diff --git a/getMesh.m b/getMesh.m index ebdc3da..be2e41b 100644 --- a/getMesh.m +++ b/getMesh.m @@ -10,6 +10,8 @@ % 'type' {'trirep'}|'triangulation' % - Specifies whether to return the mesh as a TriRep object or as a % Triangulation object +% 'limitToTriangulation' {'false'}|true +% - Specifies whether to repack the triangulation % % GETMESH Returns a face/vertex representation of the anatomical model. % Supported data types include istances of the Matlab objects Trirep and @@ -31,11 +33,14 @@ nStandardArgs = 1; % UPDATE VALUE type = 'trirep'; +limitToTriangulation = false; if nargin > nStandardArgs for i = 1:2:nargin-nStandardArgs - switch varargin{i} + switch lower(varargin{i}) case 'type' type = varargin{i+1}; + case 'limittotriangulation' + limitToTriangulation - varargin{i+1}; end end end @@ -43,8 +48,15 @@ error(['OPENEP/GETMESH: Value: ' type ' for parameter: type not recognised']); end -FV.vert = userdata.surface.triRep.X; -FV.faces = userdata.surface.triRep.Triangulation; +if isa(userdata.surface.triRep, 'TriRep') + FV.vert = userdata.surface.triRep.X; + FV.faces = userdata.surface.triRep.Triangulation; +elseif isa(userdata.surface.triRep, 'triangulation') + FV.vert = userdata.surface.triRep.Points; + FV.faces = userdata.surface.triRep.ConnectivityList; +else + error('OPENEP/getMesh: userdata.surface.TriRep should be a TriRep or a triangulation object') +end switch lower(type) case 'trirep' @@ -65,4 +77,8 @@ end end +if limitToTriangulation + tr = repack(tr); +end + end diff --git a/getNormals.m b/getNormals.m new file mode 100644 index 0000000..e2f73a0 --- /dev/null +++ b/getNormals.m @@ -0,0 +1,48 @@ +function [normals, userdata] = getNormals( userdata ) +% GETNORMALS computes or returns surface normals on OpenEP anatomy +% +% Usage: +% normals = computeNormals( userdata ) +% [normals, userdata] = computeNormals( userdata ) +% Where: +% userdata - see importcarto_mem +% normals - the surface normals at each vertex +% +% COMPUTENORMALS Returns the normals but also returns a new userdata +% structure with normals saved +% +% Author: Steven Williams (2021) (Copyright) +% SPDX-License-Identifier: Apache-2.0 +% +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% faces = getFaces( userdata ) +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +if isfield(userdata.surface, 'normals') + if isempty(userdata.surface.normals) + V = getVertices(userdata); + T = getFaces(userdata); + normals = local_computeNormals(V,T); + else + normals = userdata.surface.normals; + end +else + V = getVertices(userdata); + T = getFaces(userdata); + normals = local_computeNormals(V,T); +end + +userdata.surface.normals = normals; + + function N = local_computeNormals(V, T) + N = compute_vertex_normals(V, T, 6, 'norm'); + end + +end diff --git a/getVertices.m b/getVertices.m index 7f1fced..22a9613 100644 --- a/getVertices.m +++ b/getVertices.m @@ -41,7 +41,13 @@ end end -vertices = userdata.surface.triRep.X; +if isa(userdata.surface.triRep, 'TriRep') + vertices = userdata.surface.triRep.X; +elseif isa(userdata.surface.triRep, 'triangulation') + vertices = userdata.surface.triRep.Points; +else + error('OPENEP/getVertices: userdata.surface.TriRep should be a TriRep or a triangulation object') +end [~, isVertUsed] = repack(getMesh(userdata)); if used diff --git a/plotTag.m b/plotTag.m index a86b772..de85467 100644 --- a/plotTag.m +++ b/plotTag.m @@ -3,9 +3,10 @@ % % Usage: % h = plotTag( userdata, varargin ) +% % Where: % userdata - see importcarto_mem -% h(i) - is an array of handles referencing the plotted surfaces +% h - an array of handles representing the plotted surface % % PLOTTAG accepts the following parameter-value pairs % 'coord' {[x y x]} @@ -17,6 +18,7 @@ % - The color of the tag to draw % 'size' % - The size of the tag to draw +% 'type' {'3d'}|'surf'|'both' % % PLOTTAG Plots tag(s) on the current map % @@ -41,6 +43,7 @@ pointnum = []; color = 'r'; r = 3; +type = '3d'; if nargin > nStandardArgs for i = 1:2:nargin-1 switch varargin{i} @@ -52,28 +55,51 @@ color = varargin{i+1}; case 'size' r = varargin{i+1}; + case 'type' + type = varargin{i+1}; end end end % TODO add error checking for input parsing +% set up the constants +SURFDIAMETER = r + 0.5; +SURFTHICKNESS = 1.5; +NUMELEM = 16; + hold on -if ~isempty(coord) - % plot coordinates - for i = 1:size(coord,1) - h(i) = plotsphere(coord(i,1), coord(i,2), coord(i,3), color, r, 16); - end +if ~isempty(pointnum) + coord = [coord; userdata.electric.egmX(pointnum,:)]; +end + +if strcmpi(type, 'surf') || strcmpi(type, 'both') + N = getNormals(userdata); + surfCoord = userdata.electric.egmSurfX(pointnum,:); + vi = findclosestvertex(getMesh(userdata, 'type', 'triangulation', 'limittotriangulation', false), surfCoord); + endPoints = surfCoord + (N(vi,:) * SURFTHICKNESS); end -if ~isempty(pointnum) - % plot pointnums - coord = userdata.electric.egmX(pointnum,:); - for i = 1:numel(pointnum) - h(i) = plotsphere(coord(i,1), coord(i,2), coord(i,3), color, r, 16); +h = []; +for i = 1:size(coord,1) + + if strcmpi(type, '3d') || strcmpi(type, 'both') + h(end+1) = plotsphere(coord(i,1), coord(i,2), coord(i,3), color, r, NUMELEM); + end + if strcmpi(type, 'surf') || strcmpi(type, 'both') + ind = length(h); + [h(ind+1) h(ind+2) h(ind+3)] = drawCylinder([surfCoord(i,1), surfCoord(i,2), surfCoord(i,3) ... + , endPoints(i,1), endPoints(i,2), endPoints(i,3) ... + , SURFDIAMETER], NUMELEM ... + , 'FaceColor', color ... + ); + set(h(2:3), 'DiffuseStrength', 0.5); end + end -set(h, 'edgecolor', 'none'); +if strcmpi(type, '3d') || strcmpi(type, 'both') + set(h, 'edgecolor', 'none'); +end end diff --git a/private/colorShell.m b/private/colorShell.m index 93bfb3c..3a1e359 100644 --- a/private/colorShell.m +++ b/private/colorShell.m @@ -250,11 +250,24 @@ function adjustSlider(~, ~, ~) magenta = [1 0 1] * .8; red = cMap(1,:); %red = [1 0 0]; + case 'cv' cMap = colormap(hot); cBarTitle = 'CV (m/s)'; magenta = cMap(end,:); red = cMap(1,:); + + case 'geodesic' + cMap = flipud(colormap(winter)); + cBarTitle = 'Distance (cm)'; + + case 'wallthickness' + cMap = flipud(jet); + cMap(1:length(cMap)/10,:) = []; + magenta = [1 0 1] * .8; + red = cMap(1,:); + cBarTitle = 'Wall Thickness (mm)'; + end if ~isempty(usrColorMap) cMap = usrColorMap; diff --git a/private/drawFreeBoundary.m b/private/drawFreeBoundary.m index 6b91c97..dc8d3f5 100644 --- a/private/drawFreeBoundary.m +++ b/private/drawFreeBoundary.m @@ -30,6 +30,14 @@ drawline = varargin{1}; end +if isa(tr, 'triangulation') + Vertices = tr.Points; +elseif isa(tr, 'TriRep') + Vertices = tr.X; +else + error('OPENEP/drawFreeBoundary.m the input "tr" should be a TriRep or triangulation object'); +end + rimEdges= freeBoundary(tr); n = size(rimEdges,1); @@ -37,9 +45,9 @@ y = nan(size(x)); z = nan(size(x)); -for i = 0:(n-1); - loc1 = tr.X(rimEdges(i+1,1),:); - loc2 = tr.X(rimEdges(i+1,2),:); +for i = 0:(n-1) + loc1 = Vertices(rimEdges(i+1,1),:); + loc2 = Vertices(rimEdges(i+1,2),:); x(3*i+1) = loc1(1); y(3*i+1) = loc1(2); From 19d6df9d642cee421643b52ef4e8b53caae29c0c Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Mon, 22 Nov 2021 10:54:40 +0000 Subject: [PATCH 15/30] Refactored distanceBetweenPoints.m and distBetweenPoints.m. This is necessary to permit use of distBetweenPoints to calculate distances independently of userdata. --- distanceBetweenPoints.m | 46 +++++------------- private/distBetweenPoints.m | 93 ++++++++++++++++++++++++++++--------- 2 files changed, 82 insertions(+), 57 deletions(-) diff --git a/distanceBetweenPoints.m b/distanceBetweenPoints.m index 5986931..c51e04a 100644 --- a/distanceBetweenPoints.m +++ b/distanceBetweenPoints.m @@ -1,11 +1,13 @@ function distance = distanceBetweenPoints(userdata, P1, P2, varargin) -% DISTANCEBETWEENPOINTS Returns the distance from A to B. +% DISTANCEBETWEENPOINTS Returns the distance between electrogram recording +% locations P1 and P2. +% % Usage: % distance = distanceBetweenPoints(userdata, A, B) % Where: % userdata - see importcarto_mem -% P1 - is the first point -% P2 - is the second point +% P1 - is the index of the first electrogram point(s) +% P2 - is the index of the second electrogram point(s) % % DISTBETWEENPOINTS accepts the following parameter-value pairs % 'method' {'linear'} | 'geodesic' @@ -56,10 +58,8 @@ % get the co-ordinates of the points A = userdata.electric.egmX(P1,:); B = userdata.electric.egmX(P2,:); - - % calculate the linear distance - diffsq = (A - B).^2; - distance = sqrt(sum(diffsq, 2)); + + distance = distBetweenPoints(A, B, 'method', 'linear'); % plot a figure if plot @@ -75,33 +75,8 @@ % get the co-ordinates of the surface points A = userdata.electric.egmSurfX(P1,:); B = userdata.electric.egmSurfX(P2,:); - - % calculate the geodesic distance (repacking and reducepatch - % necessary to prevent ?memory problem in exact geodesic) - V = userdata.surface.triRep.X; - F = userdata.surface.triRep.Triangulation; - [faces,vertices] = reducepatch(F,V,size(F,1)); - - % now find the closest vertex index to A and B - P1new = findclosestvertex(vertices, A); - P2new = findclosestvertex(vertices, B); - - % geodesic algorithm - mesh = geodesic_new_mesh(vertices,faces); - algorithm = geodesic_new_algorithm(mesh, 'exact'); - source_point = {geodesic_create_surface_point('vertex',P1new,vertices(P1new,:))}; - geodesic_propagate(algorithm, source_point); % the most time-consuming step - - % find a shortest path from source to target - destination = geodesic_create_surface_point('vertex',P2new,vertices(P2new,:)); - path = geodesic_trace_back(algorithm, destination); - - % find distances - [x,y,z] = extract_coordinates_from_path(path); - distance = sum(sqrt(diff(x).^2 + diff(y).^2 + diff(z).^2)); - - %delete all meshes and algorithms - geodesic_delete; + + [distance, pathCoordinates] = distBetweenPoints(A, B, 'method', 'geodesic', 'userdata', userdata); if plot hSurf = drawMap(userdata, 'type', 'none', 'orientation', 'pa'); @@ -109,6 +84,9 @@ hold on plotTag(userdata, 'coord', A); plotTag(userdata, 'coord', B); + x = pathCoordinates(:,1); + y = pathCoordinates(:,2); + z = pathCoordinates(:,3); plot3(x*1.001,y*1.001,z*1.001,'k-','LineWidth',2); %plot path end end diff --git a/private/distBetweenPoints.m b/private/distBetweenPoints.m index 54b1acb..c92fcfe 100644 --- a/private/distBetweenPoints.m +++ b/private/distBetweenPoints.m @@ -1,24 +1,19 @@ -function D = distBetweenPoints(A, B, varargin) -% ********************************************************************** -% * The contents of this package are copyright (c) King's College London -% * -% * They may not be reproduced, distributed, modified or sold for any -% * purpose -% * -% * Author: Dr S. E. Williams -% * Address: Division of Imaging Sciences and Biomedical Engineering, -% * King's College London -% * St Thomas' Hospital -% * -% * March 2014 -% ********************************************************************** -% +function [D, pathCoordinates] = distBetweenPoints(A, B, varargin) % DISTBETWEENPOINTS Returns the distance from A to B. % Usage: % D = DISTBETWEENPOINTS(A, B) % Where: % A - is the first point(s) % B - is the second point(s) +% D - distance between the points +% pathCoordinates - co-ordinates of the geodesic path, empty if 'method' +% is 'linear' +% +% DISTBETWEENPOINTS accepts the folloiwng parameter-value pairs +% 'method' {'linear'} | 'geodesic' +% - Specifies whether to calcualte linear or geodesic distances +% 'userdata' {[]} | userdata +% - Specifies the OpenEP userdata structure for geodesic calculations % % DISTBETWEENPOINTS returns the distance from A to B. A and B are specified % as row vectors [x, y, z] or matrices, with rows representing different @@ -37,15 +32,32 @@ % code % --------------------------------------------------------------- -if nargin==3 - warn = varargin{1}; -else - warn = false; +% set up global variables +global geodesic_library; +geodesic_library = 'geodesic_matlab_api'; + +% parse input arguments +nStandardArgs = 3; +method = 'linear'; +warn = false; +userdata = []; +if nargin > nStandardArgs + for i = 1:2:nargin-nStandardArgs + switch varargin{i} + case 'method' + method = varargin{i+1}; + case 'verbose' + warn = varargin{i+1}; + case 'userdata' + userdata = varargin{i+1}; + end + end end -if size(A,2) ~= size(B,2); + +if size(A,2) ~= size(B,2) error('DISTBETWEENPOINTS: dimensions of points in A and B must be equal'); end -if size(A,1) ~= size(B,1); +if size(A,1) ~= size(B,1) if size(A,1) ~= 1 error('DISTBETWEENPOINTS: if npoints in A and B are different A must specify one and only one point'); else @@ -56,8 +68,43 @@ end end -diffsq = (A - B).^2; +switch lower(method) + case 'linear' + diffsq = (A - B).^2; + D = sqrt(sum(diffsq, 2)); + pathCoordinates = []; + + case 'geodesic' + if isempty(userdata) + error('OPENEP/DISTBETWEENPOINTS: you must specify userdata for geodesic distance calculations') + end + % calculate the geodesic distance (repacking and reducepatch + % necessary to prevent ?memory problem in exact geodesic) + V = userdata.surface.triRep.X; + F = userdata.surface.triRep.Triangulation; + [faces,vertices] = reducepatch(F,V,size(F,1)); + + % now find the closest vertex index to A and B + P1new = findclosestvertex(vertices, A); + P2new = findclosestvertex(vertices, B); + + % geodesic algorithm + mesh = geodesic_new_mesh(vertices,faces); + algorithm = geodesic_new_algorithm(mesh, 'exact'); + source_point = {geodesic_create_surface_point('vertex',P1new,vertices(P1new,:))}; + geodesic_propagate(algorithm, source_point); % the most time-consuming step + + % find a shortest path from source to target + destination = geodesic_create_surface_point('vertex',P2new,vertices(P2new,:)); + path = geodesic_trace_back(algorithm, destination); + + % find distances + [x,y,z] = extract_coordinates_from_path(path); + D = sum(sqrt(diff(x).^2 + diff(y).^2 + diff(z).^2)); -D = sqrt(sum(diffsq, 2)); + pathCoordinates = [x y z]; + + %delete all meshes and algorithms + geodesic_delete; end From 7927b0d35e5ab849ee5a8dff1faca9c685f95984 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Mon, 22 Nov 2021 12:09:17 +0000 Subject: [PATCH 16/30] added uncertaintyZones.m function --- uncertaintyZones.m | 105 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 uncertaintyZones.m diff --git a/uncertaintyZones.m b/uncertaintyZones.m new file mode 100644 index 0000000..91726bd --- /dev/null +++ b/uncertaintyZones.m @@ -0,0 +1,105 @@ +function [ points, weights ] = uncertaintyZones(userdata, loc, regionDefinition, varargin) +% UNCERTAINTYZONES identifies points of interest around locations of +% interest +% +% Usage: +% [ points, weights ] = uncertaintyZones( userdata, loc, regionDefinition ) +% Where: +% userdata - see importcarto_mem +% loc - the co-ordinates of the points we are interested in, n x 3 +% regionDefinition - structure with the following fields: +% .shape - the shape, only sphere is implemented +% .params - parameters describing the shape, only radius is applicable +% points - logical array of size n x m +% weights - numerical array of size n x m +% +% UNCERTAINTYZONES accepts the following parameter-value pairs +% 'type' {'surface'}|'electrogram' +% 'weightalgorithm' {'linear'} +% +% UNCERTAINTYZONES is used to identify points of interest around locations +% of interest. For example, UNCERTAINTYZONES can be used to access data +% such as electrogram voltage data within a region around each mapping +% point. +% +% The function works by calculating a radius aroudn each location of +% interest (loc), and identifying every data point that falls within that +% sphere. A weighting is calculated for each data point of interest based +% on the weighting algorithm, currently the only option is _linear_. +% +% UNCERTAINTYZONES can work on surface data or electrogram data. In the +% case of surface data, the mesh nodes are taken as the set of all possible +% data locations. In the case of electrogram data, the electrogram +% recordings locations are taken as the set of all possible data locations. +% +% The results are returned as two arrays, points and weights. The size of +% both is n x m, where n is the number of locations queried, and m is the +% number of all data locations available. Points takes values of 0 at all +% data locations outwith the scope of the region definition and 1 at all +% locations within the scope of the region definition. Weights has a +% numberical value representing the weights at all locations, in the range +% [0 1]. By definition, any data locations outwith the region definition +% will have a weight of zero. +% +% Author: Steven Williams (2021) +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% test code +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% Parse inputs +nStandardArgs = 3; % UPDATE VALUE +type = 'map'; +weightAlgorithm = 'linear'; +if nargin > nStandardArgs + for i = 1:2:nargin-nStandardArgs + switch varargin{i} + case 'type' + type = varargin{i+1}; + case 'weightalgorithm' + weightAlgorithm = varargin{i+1}; + end + end +end + +% get the potential data positions +switch lower(type) + case 'map' + X = getVertices(userdata); % todo - how do we handle points not referenced by the triangulation + case 'electrogram' + X = getElectrogramX(userdata); +end + +% identify the points of interest +m = size(X,1); +n = size(loc,1); +points = NaN(n, m); +weights = zeros(n,m); + +switch lower(regionDefinition.shape) + case 'sphere' + radius = regionDefinition.params(1); + for n = 1:size(loc,1) + distances = distBetweenPoints(loc(n,:),X,'method','linear'); + points(n,distances0) = 0; + weights(n,:) = abs(temp); + end +end + +end + + + + + + + From e3e11160451a0e70242d8f170cd67b2a1c76b112 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Mon, 22 Nov 2021 16:57:15 +0000 Subject: [PATCH 17/30] Added a function to downsample an OpenEP dataset --- downSampleDataset.m | 71 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 downSampleDataset.m diff --git a/downSampleDataset.m b/downSampleDataset.m new file mode 100644 index 0000000..97c8cec --- /dev/null +++ b/downSampleDataset.m @@ -0,0 +1,71 @@ +function [ userdata_ds ] = downSampleDataset( userdata, sampleDensity ) +% DOWNSAMPLEDATASET Returns a new OpenEP dataset at the required sampling +% density +% +% Usage: +% [ userdata_ds ] = DOWNSAMPLEDATASET( userdata, sampleDensity ) +% Where: +% userdata - the input +% userdata_ds - the output +% +% DOWNSAMPLEDATASET Does not check the distribution of points but just +% returns the required number of points per map +% +% Author: Steven Williams (2021) +% Modifications - +% +% Info on Code Testing: +% --------------------------------------------------------------- +% test code +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% Work out the surface area and original density +surfaceArea = getArea(userdata); +numPts = getNumPts(userdata); +densityOriginal = numPts/surfaceArea; +disp(['Original density = ' num2str(densityOriginal) 'points/cm2']) + +% Calculate and identify the required number and points +requiredNumberOfPoints = ceil(surfaceArea * sampleDensity); +if requiredNumberOfPoints > numPts + error('OPENEP/DOWNSAMPLEDATASET: Original sampling density is insufficient to fulfill this request'); +end +indices = randperm(numPts); +iKeep = indices(1:requiredNumberOfPoints); +iRemove = indices(requiredNumberOfPoints+1:end); + +% Remove the data we do not need +userdata_ds = userdata; +userdata_ds.electric.tags(iRemove,:) = []; +userdata_ds.electric.names(iRemove,:) = []; +userdata_ds.electric.electrodeNames_bip(iRemove,:) = []; +userdata_ds.electric.egmX(iRemove,:) = []; +userdata_ds.electric.egm(iRemove,:) = []; +userdata_ds.electric.electrodeNames_uni(iRemove,:) = []; +userdata_ds.electric.egmUniX(iRemove,:,:) = []; +userdata_ds.electric.egmUni(iRemove,:,:) = []; +userdata_ds.electric.egmRef(iRemove,:) = []; +userdata_ds.electric.ecg(iRemove,:) = []; +userdata_ds.electric.annotations.woi(iRemove,:) = []; +userdata_ds.electric.annotations.referenceAnnot(iRemove,:) = []; +userdata_ds.electric.annotations.mapAnnot(iRemove,:) = []; +userdata_ds.electric.voltages.bipolar(iRemove,:) = []; +userdata_ds.electric.voltages.unipolar(iRemove,:) = []; +userdata_ds.electric.impedances.time(:,iRemove) = []; +userdata_ds.electric.impedances.value(:,iRemove) = []; +userdata_ds.electric.egmSurfX(iRemove,:) = []; +userdata_ds.electric.barDirection(iRemove,:) = []; + +% Work out the new density +numPts_ds = getNumPts(userdata_ds); +density_ds = numPts_ds/surfaceArea; +disp(['New density = ' num2str(density_ds) 'points/cm2']) + +% Add a note +userdata_ds.notes{end+1} = [date ': downsampled'] + +end \ No newline at end of file From 23b444564f23ce1388cca69ca0bfa218d0b1bc58 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 23 Nov 2021 06:46:53 +0000 Subject: [PATCH 18/30] small change to output --- uncertaintyZones.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uncertaintyZones.m b/uncertaintyZones.m index 91726bd..612aa6a 100644 --- a/uncertaintyZones.m +++ b/uncertaintyZones.m @@ -79,7 +79,7 @@ % identify the points of interest m = size(X,1); n = size(loc,1); -points = NaN(n, m); +points = zeros(n, m); weights = zeros(n,m); switch lower(regionDefinition.shape) From 6731d6868cd4468e6cb3760fedb47f1785cb4f69 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 23 Nov 2021 08:37:10 +0000 Subject: [PATCH 19/30] Added functionality for rendering weights into drawMap --- drawMap.m | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/drawMap.m b/drawMap.m index 5c88d42..4b865d8 100644 --- a/drawMap.m +++ b/drawMap.m @@ -74,7 +74,7 @@ orientation = inputs.orientation; DISTANCETHRESH = inputs.colorfillthreshold; -if ~any(strcmpi(type, {'act', 'bip', 'force', 'uni', 'none', 'cv', 'geodesic', 'wallthickness'})) +if ~any(strcmpi(type, {'act', 'bip', 'force', 'uni', 'none', 'cv', 'geodesic', 'wallthickness', 'weights'})) error('DRAWMAP: Invalid parameter-value pair; type must be one of act, bip, force, uni, none, cv or geodesic') end @@ -144,6 +144,15 @@ else error('OPENEP/drawMap: automatic calculation of wall thickness data not supported') end + case 'weights' + DATATYPE = 'weights'; + if ~isempty(DATAMANUAL) + DATA = DATAMANUAL; + else + [~, weights] = uncertaintyZones(userdata, getElectrogramX(userdata), []); + DATA = max(weights, [], 2); + end + end if ~strcmpi(type, 'none') From cb30819aea2aeab1c5a302ede7280c134d8dd85a Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 23 Nov 2021 08:37:32 +0000 Subject: [PATCH 20/30] Minor update to documentation fro getElectrogramX --- getElectrogramX.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/getElectrogramX.m b/getElectrogramX.m index 3dd9921..9a1ad99 100644 --- a/getElectrogramX.m +++ b/getElectrogramX.m @@ -2,7 +2,7 @@ % GETELECTROGRAMX Returns the the electrode recording positions % % Usage: -% C = getCentreOfMass( userdata, varargin ) +% [X, surfX] = getElectrogramX( userdata, varargin ) % Where: % userdata - see importcarto_mem % X - the 3D Cartesian co-ordinates From aba7636db9c4efa4cf53308003dd56bb1e7e3bfa Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 23 Nov 2021 08:38:00 +0000 Subject: [PATCH 21/30] Change transposition of output and added a default region definition. --- uncertaintyZones.m | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/uncertaintyZones.m b/uncertaintyZones.m index 612aa6a..98e1583 100644 --- a/uncertaintyZones.m +++ b/uncertaintyZones.m @@ -68,6 +68,12 @@ end end +% Create default regionDefinition +if isempty(regionDefinition) + regionDefinition.shape = 'sphere'; + regionDefinition.params = 5; +end + % get the potential data positions switch lower(type) case 'map' @@ -95,6 +101,9 @@ end end +points = points'; +weights = weights'; + end From 7b2293596b87f5a42a97745640d6bda139fa3c36 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 23 Nov 2021 08:38:15 +0000 Subject: [PATCH 22/30] adding weights --- private/colorShell.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/private/colorShell.m b/private/colorShell.m index 3a1e359..348f033 100644 --- a/private/colorShell.m +++ b/private/colorShell.m @@ -267,6 +267,10 @@ function adjustSlider(~, ~, ~) magenta = [1 0 1] * .8; red = cMap(1,:); cBarTitle = 'Wall Thickness (mm)'; + + case 'weights' + cMap = colormap(parula); + cBarTitle = 'Normalised distance weights'; end if ~isempty(usrColorMap) From d044fc60443c23a6d77e5f95e657b80f1dadc970 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 23 Nov 2021 08:38:55 +0000 Subject: [PATCH 23/30] Only compute normals if they do not already exist. --- plotTag.m | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/plotTag.m b/plotTag.m index de85467..3cb6f47 100644 --- a/plotTag.m +++ b/plotTag.m @@ -1,4 +1,4 @@ -function h = plotTag( userdata, varargin ) +function [h vi] = plotTag( userdata, varargin ) % PLOTTAG Plots tag(s) on the current map % % Usage: @@ -7,6 +7,7 @@ % Where: % userdata - see importcarto_mem % h - an array of handles representing the plotted surface +% vi - indices of the surface projection if 'type' is 'surf' or 'both' % % PLOTTAG accepts the following parameter-value pairs % 'coord' {[x y x]} @@ -72,10 +73,21 @@ coord = [coord; userdata.electric.egmX(pointnum,:)]; end +vi = []; if strcmpi(type, 'surf') || strcmpi(type, 'both') - N = getNormals(userdata); - surfCoord = userdata.electric.egmSurfX(pointnum,:); - vi = findclosestvertex(getMesh(userdata, 'type', 'triangulation', 'limittotriangulation', false), surfCoord); + if ~isfield(userdata.surface, 'normals') + N = getNormals(userdata); + else + N = userdata.surface.normals; + end + if ~isempty(pointnum) + surfCoord = userdata.electric.egmSurfX(pointnum,:); + vi = findclosestvertex(getMesh(userdata, 'type', 'triangulation', 'limittotriangulation', false), surfCoord); + else + mesh = getMesh(userdata, 'type', 'triangulation', 'limittotriangulation', false); + vi = findclosestvertex(mesh, coord); + surfCoord = mesh.Points(vi,:); + end endPoints = surfCoord + (N(vi,:) * SURFTHICKNESS); end From 0ae8acb4f1d9fb54689a56ab3d9a0511634c01c8 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Thu, 16 Dec 2021 14:42:00 +0000 Subject: [PATCH 24/30] corrected problem with gettMesh.m arising from the previous merge --- getMesh.m | 45 ++++++++++++++------------------------------- 1 file changed, 14 insertions(+), 31 deletions(-) diff --git a/getMesh.m b/getMesh.m index e2ff7b7..44835f0 100644 --- a/getMesh.m +++ b/getMesh.m @@ -7,18 +7,14 @@ % tr - a TriRep, or Triangulation, object % % GETMESH accepts the following parameter-value pairs -% 'type' {'trirep'}|'triangulation' +% 'type' {'trirep'}|'triangulation'|'struct' % - Specifies whether to return the mesh as a TriRep object or as a -% Triangulation object -<<<<<<< HEAD +% Triangulation object or as a struct % 'limitToTriangulation' {'false'}|true % - Specifies whether to repack the triangulation -======= -% 'repack' {false}|true ->>>>>>> origin/develop % -% GETMESH Returns a face/vertex representation of the anatomical model. -% Supported data types include istances of the Matlab objects Trirep and +% GETMESH Returns a face/vertex representation of the anatomical model. +% Supported data types include istances of the Matlab objects Trirep and % Triangulation. % % Author: Steven Williams (2020) (Copyright) @@ -37,45 +33,33 @@ nStandardArgs = 1; % UPDATE VALUE type = 'trirep'; -<<<<<<< HEAD limitToTriangulation = false; -======= dorepack = false; ->>>>>>> origin/develop + if nargin > nStandardArgs for i = 1:2:nargin-nStandardArgs switch lower(varargin{i}) case 'type' type = varargin{i+1}; -<<<<<<< HEAD case 'limittotriangulation' - limitToTriangulation - varargin{i+1}; -======= + limitToTriangulation = varargin{i+1}; case 'repack' dorepack = varargin{i+1}; ->>>>>>> origin/develop end end end -if ~any(strcmpi({'trirep' 'triangulation'}, type)) +if ~any(strcmpi({'trirep' 'triangulation' 'struct'}, type)) error(['OPENEP/GETMESH: Value: ' type ' for parameter: type not recognised']); end -<<<<<<< HEAD -======= % Check the type of surface data stored in the OpenEP datastructure and % accordingly access the triangulation and the points ->>>>>>> origin/develop if isa(userdata.surface.triRep, 'TriRep') FV.vert = userdata.surface.triRep.X; FV.faces = userdata.surface.triRep.Triangulation; elseif isa(userdata.surface.triRep, 'triangulation') FV.vert = userdata.surface.triRep.Points; FV.faces = userdata.surface.triRep.ConnectivityList; -<<<<<<< HEAD -else - error('OPENEP/getMesh: userdata.surface.TriRep should be a TriRep or a triangulation object') -======= elseif isa(userdata.surface.triRep, 'struct') if ~isfield(userdata.surface.triRep, 'X') error('OPENEP/GETMESH: invalid data. userdata.surface.triRep must be one of: TriRep, triangulation or struct with fields .X and .Triangulation'); @@ -85,7 +69,6 @@ end FV.vert = userdata.surface.triRep.X; FV.faces = userdata.surface.triRep.Triangulation; ->>>>>>> origin/develop end switch lower(type) @@ -101,15 +84,15 @@ else tr = triangulation(FV.faces, FV.vert(:,1), FV.vert(:,2), FV.vert(:,3)); end + case 'struct' + if isa(userdata.surface.triRep, 'struct') + tr = userdata.surface.triRep; + else + tr.X = FV.vert; + tr.Points = FV.faces; + end end -<<<<<<< HEAD if limitToTriangulation tr = repack(tr); end - -======= -if dorepack - tr = repack(tr); ->>>>>>> origin/develop -end From 6193815350ecb932686756fe56d3240dd3b94027 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Thu, 16 Dec 2021 14:43:00 +0000 Subject: [PATCH 25/30] corrected propblem with drawFreeBoundary.m arising from previous merge --- private/drawFreeBoundary.m | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/private/drawFreeBoundary.m b/private/drawFreeBoundary.m index 48af03f..92f4a1b 100644 --- a/private/drawFreeBoundary.m +++ b/private/drawFreeBoundary.m @@ -31,40 +31,30 @@ end if isa(tr, 'triangulation') - Vertices = tr.Points; + X = tr.Points; elseif isa(tr, 'TriRep') - Vertices = tr.X; + X = tr.X; +elseif isa(tr, 'struct') + X = tr.X; else - error('OPENEP/drawFreeBoundary.m the input "tr" should be a TriRep or triangulation object'); + error('OPENEP/drawFreeBoundary.m the input "tr" should be a TriRep object, a triangulation object or a structure'); end -rimEdges= freeBoundary(tr); +rimEdges = freeBoundary(tr); n = size(rimEdges,1); x = nan(3*n,1); y = nan(size(x)); z = nan(size(x)); -<<<<<<< HEAD -for i = 0:(n-1) - loc1 = Vertices(rimEdges(i+1,1),:); - loc2 = Vertices(rimEdges(i+1,2),:); -======= -if isa(tr, 'TriRep') - X = tr.X; -elseif isa(tr, 'triangulation') - X = tr.Points; -end - for i = 0:(n-1) loc1 = X(rimEdges(i+1,1),:); loc2 = X(rimEdges(i+1,2),:); ->>>>>>> origin/develop - + x(3*i+1) = loc1(1); y(3*i+1) = loc1(2); z(3*i+1) = loc1(3); - + x(3*i+2) = loc2(1); y(3*i+2) = loc2(2); z(3*i+2) = loc2(3); @@ -86,4 +76,4 @@ h.Z = z; end -end \ No newline at end of file +end From 0379fa70ea921015b336e919cff0c9531520b248 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Tue, 11 Jan 2022 17:14:13 +0000 Subject: [PATCH 26/30] comment added --- uncertaintyZones.m | 1 + 1 file changed, 1 insertion(+) diff --git a/uncertaintyZones.m b/uncertaintyZones.m index 98e1583..ccca2e4 100644 --- a/uncertaintyZones.m +++ b/uncertaintyZones.m @@ -72,6 +72,7 @@ if isempty(regionDefinition) regionDefinition.shape = 'sphere'; regionDefinition.params = 5; + disp('Default region definition used'); end % get the potential data positions From 2baaf838af668aa342171f3ea090ca4be3046d0c Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Wed, 10 May 2023 06:20:33 +0100 Subject: [PATCH 27/30] updated changelog for uncertaintyzones --- changelog.md | 1 + 1 file changed, 1 insertion(+) diff --git a/changelog.md b/changelog.md index 75bc6ac..ef4c07c 100644 --- a/changelog.md +++ b/changelog.md @@ -19,6 +19,7 @@ This section documents changes which have been merged into develop and will form naming conventions of catheters in Carto3! - comparestructure.m allows two structures to be compared, eg userdata from two different versions of OpenEP +- uncertaintyzones.m calculates regional analysis around mapping points ### Changed - Surface models stored as structures rather than TriRep objects for compatibility with OpenEP From b8794059ebf625b2bb897fee4be03a6f9d3071af Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Wed, 10 May 2023 06:36:42 +0100 Subject: [PATCH 28/30] added condition for structure to getVertices --- getVertices.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/getVertices.m b/getVertices.m index 22a9613..b769570 100644 --- a/getVertices.m +++ b/getVertices.m @@ -45,8 +45,10 @@ vertices = userdata.surface.triRep.X; elseif isa(userdata.surface.triRep, 'triangulation') vertices = userdata.surface.triRep.Points; +elseif isa(userdata.surface.triRep, 'struct') + vertices = userdata.surface.triRep.X; else - error('OPENEP/getVertices: userdata.surface.TriRep should be a TriRep or a triangulation object') + error('OPENEP/getVertices: userdata.surface.TriRep should be a Structure, TriRep or a triangulation object') end [~, isVertUsed] = repack(getMesh(userdata)); From 5f6a5af91f4b52660d15665c52d2e8cc55d54ec0 Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Wed, 10 May 2023 06:40:43 +0100 Subject: [PATCH 29/30] added some further documentation for uncertaintyzones.m --- uncertaintyZones.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/uncertaintyZones.m b/uncertaintyZones.m index ccca2e4..9fbad77 100644 --- a/uncertaintyZones.m +++ b/uncertaintyZones.m @@ -46,7 +46,11 @@ % % Info on Code Testing: % --------------------------------------------------------------- -% test code +% [~, surfX] = getElectrogramX(userdata); +% regionDefinition.shape = 'sphere'; +% regionDefinition.params = 6; +% [points, weights] = uncertaintyZones(userdata, surfX(1,:), regionDefinition); +% drawMap(userdata, 'data', points) % --------------------------------------------------------------- % % --------------------------------------------------------------- From bc18ccc9dc2e45beeddaa009fb09524a354a177f Mon Sep 17 00:00:00 2001 From: Steven Williams Date: Mon, 28 Aug 2023 14:35:57 +0100 Subject: [PATCH 30/30] minor changes --- batchConvert.m | 36 +++++++++++++++++++++++++++++++++--- batchProcess.m | 2 ++ getFaces.m | 4 +++- 3 files changed, 38 insertions(+), 4 deletions(-) diff --git a/batchConvert.m b/batchConvert.m index d96b28a..90364c7 100644 --- a/batchConvert.m +++ b/batchConvert.m @@ -1,9 +1,10 @@ -function batchConvert(inputDir, outputDir) +function batchConvert(inputDir, outputDir, varargin) % batchConvert Converts original style OpenEP datasets into the new format. % % In the newer format, userdata.surface.triRep is a regular struct rather % than a TriRep. This is to allow for interopability with OpenEP-Py, which -% cannot load TriRep objects. +% cannot load TriRep objects. This function also allows options to +% pre-compute normals. % % Usage: % batchConvert(inputDir, outputDir) @@ -15,11 +16,32 @@ function batchConvert(inputDir, outputDir) % will be the same as those in inputDir. % % Author: Paul Smith (2021) (Copyright) +% SPDX-License-Identifier: Apache-2.0 +% % Modifications - % Steven Williams (2022) Updated documentation and added notes +% Steven Williams (2023) Added option to pre-compute normals % -% SPDX-License-Identifier: Apache-2.0 +% Info on Code Testing: +% --------------------------------------------------------------- +% batchConvert(inputDir, outputDir, 'computenormals', true) +% --------------------------------------------------------------- % +% --------------------------------------------------------------- +% code +% --------------------------------------------------------------- + +% parse input arguments +nStandardArgs = 2; +computeNormals = false; +if nargin > nStandardArgs + for i = nStandardArgs:2:nargin-1-nStandardArgs + switch lower(varargin{i}) + case 'computenormals' + computeNormals = varargin{i+1}; + end + end +end % Create the output folder if it does not exist if ~exist(outputDir, 'dir') @@ -44,6 +66,14 @@ function batchConvert(inputDir, outputDir) % Store comment about what we have done userdata.notes{end+1} = [date ': data set converted using batchConvert.m']; + if computeNormals + disp('precomputing normals') + [~, userdata] = getNormals(userdata); + + % Store comment about what we have done + userdata.notes{end+1} = [date ': normals added using batchConvert.m']; + end + % We save as -v7 because it's faster to load in OpenEP-py than -v7.3, % and the saved file is significantly smaller compared to -v6 files. outputFile = [outputDir filesep() allFiles{i}]; diff --git a/batchProcess.m b/batchProcess.m index 017a546..53e9eb4 100644 --- a/batchProcess.m +++ b/batchProcess.m @@ -19,6 +19,8 @@ % Iterate through each file and perform some OpenEP functions for i = 1:numel(allFiles) + + disp(['working on file ... ' allFiles{i}]) load([working_dir filesep() allFiles{i}]); diff --git a/getFaces.m b/getFaces.m index 65775a5..1065149 100644 --- a/getFaces.m +++ b/getFaces.m @@ -27,8 +27,10 @@ faces = userdata.surface.triRep.Triangulation; elseif isa(userdata.surface.triRep, 'triangulation') faces = userdata.surface.triRep.ConnectivityList; +elseif isa(userdata.surface.triRep, 'struct') + faces = userdata.surface.triRep.Triangulation; else - error('OPENEP/getFaces: userdata.surface.TriRep should be a TriRep or a triangulation object') + error('OPENEP/getFaces: userdata.surface.TriRep should be a TriRep, a triangulation object or a structure') end end