diff --git a/UML/interpolation architecture.drawio b/UML/interpolation architecture.drawio new file mode 100644 index 0000000..bfa6d84 --- /dev/null +++ b/UML/interpolation architecture.drawio @@ -0,0 +1,114 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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/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 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_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..21d3f14 --- /dev/null +++ b/doCvMapping_Gradient.m @@ -0,0 +1,75 @@ +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---------------------------- +% 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 +% --------------------------------------------------------------- + +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 conduction velocities +sgrad = grad .* grad; +dp = sum(sgrad,2); +cv = 1./sqrt(dp); + +% 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']); + +% 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']); + +% interpolate +interpCv = int.interpolate(cvX, cv, getVertices(userdata)); + +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..1a8079e --- /dev/null +++ b/doCvMapping_RadialBasis.m @@ -0,0 +1,54 @@ +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 +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'; + +% 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/drawMap.m b/drawMap.m index 51364e1..6f886fa 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', 'weights'})) + error('DRAWMAP: Invalid parameter-value pair; type must be one of act, bip, force, uni, none, cv or geodesic') end @@ -124,6 +124,35 @@ 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 + 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') @@ -148,7 +177,7 @@ ); else disp('no data to color the shell with') - end + end end % Adjust the light/material diff --git a/getConductionVelocity.m b/getConductionVelocity.m index a053e14..ca1e419 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'`; 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 +% 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,49 @@ % code % --------------------------------------------------------------- -lats = userdata.electric.annotations.mapAnnot(getMappingPointsWithinWoI(userdata)); -X = userdata.electric.egmX(getMappingPointsWithinWoI(userdata),:); -[~,~,~,~,cvdata] = RBFConductionVelocity(lats', X', getMesh(userdata).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 = openEpDataInterpolator(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, int); + + case 'radialbasis' + [cv, cvX, interpCv] = doCvMapping_RadialBasis(userdata, int); + + 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/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/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 diff --git a/getFaces.m b/getFaces.m index e743347..1065149 100644 --- a/getFaces.m +++ b/getFaces.m @@ -23,6 +23,14 @@ % code % --------------------------------------------------------------- -faces = getMesh(userdata, 'type', 'triangulation').ConnectivityList; +if isa(userdata.surface.triRep, 'TriRep') + 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, a triangulation object or a structure') +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 55535cf..b769570 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,8 +30,29 @@ % code % --------------------------------------------------------------- -tr = getMesh(userdata, 'type', 'triangulation'); -vertices = tr.Points; -[~, isVertUsed] = repack(tr); +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 +if isa(userdata.surface.triRep, 'TriRep') + 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 Structure, TriRep or a triangulation object') end +[~, isVertUsed] = repack(getMesh(userdata)); + +if used + vertices = vertices(isVertUsed,:); + isVertUsed = true(size(vertices)); +end \ No newline at end of file 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/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 new file mode 100644 index 0000000..7c09565 --- /dev/null +++ b/openEpDataInterpolator.m @@ -0,0 +1,132 @@ +classdef openEpDataInterpolator + % 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 + distanceThreshold + end + + properties (Access = private) + call_interpolator + end + + methods + % Contructor + function obj = openEpDataInterpolator(method) + + possibleMethods = {'scatteredinterpolant' ... + , 'localsmoothing' ... + , 'radialbasis' ... + }; + + defaultMethod = 'scatteredinterpolant'; + + 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 + end + end + end + + end + + function f_q = interpolate(obj, x, f_x, q) + + f_q = obj.call_interpolator.interpolate(x, f_x, q); + + 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/plotTag.m b/plotTag.m index a86b772..3cb6f47 100644 --- a/plotTag.m +++ b/plotTag.m @@ -1,11 +1,13 @@ -function h = plotTag( userdata, varargin ) +function [h vi] = plotTag( userdata, varargin ) % PLOTTAG Plots tag(s) on the current map % % 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 +% vi - indices of the surface projection if 'type' is 'surf' or 'both' % % PLOTTAG accepts the following parameter-value pairs % 'coord' {[x y x]} @@ -17,6 +19,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 +44,7 @@ pointnum = []; color = 'r'; r = 3; +type = '3d'; if nargin > nStandardArgs for i = 1:2:nargin-1 switch varargin{i} @@ -52,28 +56,62 @@ 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); +if ~isempty(pointnum) + coord = [coord; userdata.electric.egmX(pointnum,:)]; +end + +vi = []; +if strcmpi(type, 'surf') || strcmpi(type, 'both') + 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 -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 3a839ea..348f033 100644 --- a/private/colorShell.m +++ b/private/colorShell.m @@ -250,18 +250,33 @@ 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 'labels' - cBarTitle = 'labels'; + + 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)'; + + case 'weights' + cMap = colormap(parula); + cBarTitle = 'Normalised distance weights'; + end if ~isempty(usrColorMap) cMap = usrColorMap; end - + % Limit the size of the colormap if length(cMap) > 256 cMap = downsample(cMap, round(length(cMap)/256)); @@ -326,7 +341,7 @@ function adjustSlider(~, ~, ~) cMapRed = repmat(red, [nBrown, 1]); figureCMap = vertcat(cMapRed, cMapMagenta); end - colormap(hAx, figureCMap); + colormap(figureCMap); end % Show the colorbar @@ -350,4 +365,4 @@ function adjustSlider(~, ~, ~) end end %adjustSlider() -end %colorShell() +end %colorShell() \ No newline at end of file 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 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 diff --git a/uncertaintyZones.m b/uncertaintyZones.m new file mode 100644 index 0000000..9fbad77 --- /dev/null +++ b/uncertaintyZones.m @@ -0,0 +1,119 @@ +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: +% --------------------------------------------------------------- +% [~, surfX] = getElectrogramX(userdata); +% regionDefinition.shape = 'sphere'; +% regionDefinition.params = 6; +% [points, weights] = uncertaintyZones(userdata, surfX(1,:), regionDefinition); +% drawMap(userdata, 'data', points) +% --------------------------------------------------------------- +% +% --------------------------------------------------------------- +% 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 + +% Create default regionDefinition +if isempty(regionDefinition) + regionDefinition.shape = 'sphere'; + regionDefinition.params = 5; + disp('Default region definition used'); +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 = zeros(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 + +points = points'; +weights = weights'; + +end + + + + + + +