Feature/uncertainty zones

UML/interpolation architecture.drawio

@@ -0,0 +1,114 @@
+<mxfile host="app.diagrams.net" modified="2021-07-02T14:53:13.040Z" agent="5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/605.1.15 (KHTML, like Gecko) Version/14.1 Safari/605.1.15" etag="9n6kZVsT_4Ia6h9Szg5Z" version="14.8.3" type="github">
+  <diagram id="jgYHGT9bPeA_kPWx9DmP" name="Page-1">
+    <mxGraphModel dx="599" dy="727" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+      <root>
+        <mxCell id="0" />
+        <mxCell id="1" parent="0" />
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-1" value="openEpDataInterpolator" style="swimlane;fontStyle=1;align=center;verticalAlign=top;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;" parent="1" vertex="1">
+          <mxGeometry x="20" y="150" width="280" height="248" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-2" value="+ method: String (public)&#xa;&#xa;+ distanceThreshold: Double (public)&#xa;&#xa;+ call_interpolator: openEpInterpolator (private)" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-1" vertex="1">
+          <mxGeometry y="26" width="280" height="94" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-3" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-1" vertex="1">
+          <mxGeometry y="120" width="280" height="8" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-4" value="+ openEpDataInterpolator(): constructor&#xa;&#xa;+ interpolate(): public&#xa;&#xa;+ call_interpolator(): private&#xa;&#xa;+ disp(): public, overridden" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-1" vertex="1">
+          <mxGeometry y="128" width="280" height="120" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-5" value="openEpScatteredInterpolant" style="swimlane;fontStyle=1;align=center;verticalAlign=top;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;" parent="1" vertex="1">
+          <mxGeometry x="420" y="480" width="190" height="124" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-6" value="+ interMethod: String&#xa;&#xa;+ exterMethod: String" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-5" vertex="1">
+          <mxGeometry y="26" width="190" height="64" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-7" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-5" vertex="1">
+          <mxGeometry y="90" width="190" height="8" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-8" value="+ interpolate(x, f_x, q): public" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-5" vertex="1">
+          <mxGeometry y="98" width="190" height="26" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-9" value="openEpLocalSmoothing" style="swimlane;fontStyle=1;align=center;verticalAlign=top;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;" parent="1" vertex="1">
+          <mxGeometry x="20" y="480" width="190" height="124" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-10" value="+ smoothingLength: Double&#xa;&#xa;+ fillWith: Double" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-9" vertex="1">
+          <mxGeometry y="26" width="190" height="64" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-11" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-9" vertex="1">
+          <mxGeometry y="90" width="190" height="8" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-12" value="+ interpolate(x, f_x, q): public" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-9" vertex="1">
+          <mxGeometry y="98" width="190" height="26" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-13" value="openEpLocalSmoothing" style="swimlane;fontStyle=1;align=center;verticalAlign=top;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;" parent="1" vertex="1">
+          <mxGeometry x="220" y="480" width="190" height="124" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-14" value="+ rbfFunction: String&#xa;&#xa;+ rfbConstant: Double" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-13" vertex="1">
+          <mxGeometry y="26" width="190" height="64" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-15" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-13" vertex="1">
+          <mxGeometry y="90" width="190" height="8" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-16" value="+ interpolate(x, f_x, q): public" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" parent="zvqzcq4YEUCvtmfPgBdI-13" vertex="1">
+          <mxGeometry y="98" width="190" height="26" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-17" value="" style="endArrow=open;endFill=1;endSize=12;html=1;exitX=1;exitY=0.5;exitDx=0;exitDy=0;entryX=0;entryY=0.25;entryDx=0;entryDy=0;" parent="1" source="zvqzcq4YEUCvtmfPgBdI-2" target="zvqzcq4YEUCvtmfPgBdI-27" edge="1">
+          <mxGeometry width="160" relative="1" as="geometry">
+            <mxPoint x="270" y="510" as="sourcePoint" />
+            <mxPoint x="442.06999999999994" y="182.01199999999994" as="targetPoint" />
+            <Array as="points" />
+          </mxGeometry>
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-22" value="Extends" style="endArrow=block;endSize=16;endFill=0;html=1;entryX=0.25;entryY=1;entryDx=0;entryDy=0;" parent="1" source="zvqzcq4YEUCvtmfPgBdI-9" target="zvqzcq4YEUCvtmfPgBdI-33" edge="1">
+          <mxGeometry y="10" width="160" relative="1" as="geometry">
+            <mxPoint x="730" y="870" as="sourcePoint" />
+            <mxPoint x="490" y="310" as="targetPoint" />
+            <mxPoint as="offset" />
+          </mxGeometry>
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-34" value="Extends" style="endArrow=block;endSize=16;endFill=0;html=1;entryX=0.5;entryY=1;entryDx=0;entryDy=0;" parent="1" source="zvqzcq4YEUCvtmfPgBdI-13" target="zvqzcq4YEUCvtmfPgBdI-33" edge="1">
+          <mxGeometry y="10" width="160" relative="1" as="geometry">
+            <mxPoint x="269.1077441077441" y="490" as="sourcePoint" />
+            <mxPoint x="550" y="320" as="targetPoint" />
+            <mxPoint as="offset" />
+          </mxGeometry>
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-35" value="Extends" style="endArrow=block;endSize=16;endFill=0;html=1;entryX=0.75;entryY=1;entryDx=0;entryDy=0;" parent="1" source="zvqzcq4YEUCvtmfPgBdI-5" target="zvqzcq4YEUCvtmfPgBdI-33" edge="1">
+          <mxGeometry y="10" width="160" relative="1" as="geometry">
+            <mxPoint x="468.08510638297867" y="490" as="sourcePoint" />
+            <mxPoint x="550" y="270" as="targetPoint" />
+            <mxPoint as="offset" />
+          </mxGeometry>
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-36" value="" style="group" parent="1" vertex="1" connectable="0">
+          <mxGeometry x="350" y="260" width="180" height="100" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-27" value="«interface»&lt;br&gt;&lt;b&gt;openEpInterpolator&lt;/b&gt;" style="html=1;" parent="zvqzcq4YEUCvtmfPgBdI-36" vertex="1">
+          <mxGeometry width="180" height="50" as="geometry" />
+        </mxCell>
+        <mxCell id="zvqzcq4YEUCvtmfPgBdI-33" value="&lt;span style=&quot;text-align: left&quot;&gt;+ interpolate(x, f_x, q): public&lt;/span&gt;" style="html=1;" parent="zvqzcq4YEUCvtmfPgBdI-36" vertex="1">
+          <mxGeometry y="50" width="180" height="50" as="geometry" />
+        </mxCell>
+        <mxCell id="4owkeiqQ8d7qHkTqZ-NT-1" value="quLATiInterpolant" style="swimlane;fontStyle=1;align=center;verticalAlign=top;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;" vertex="1" parent="1">
+          <mxGeometry x="620" y="480" width="190" height="124" as="geometry" />
+        </mxCell>
+        <mxCell id="4owkeiqQ8d7qHkTqZ-NT-2" value="+&#xa;&#xa;+" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" vertex="1" parent="4owkeiqQ8d7qHkTqZ-NT-1">
+          <mxGeometry y="26" width="190" height="64" as="geometry" />
+        </mxCell>
+        <mxCell id="4owkeiqQ8d7qHkTqZ-NT-3" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;" vertex="1" parent="4owkeiqQ8d7qHkTqZ-NT-1">
+          <mxGeometry y="90" width="190" height="8" as="geometry" />
+        </mxCell>
+        <mxCell id="4owkeiqQ8d7qHkTqZ-NT-4" value="+ interpolate(x, f_x, q): public" style="text;strokeColor=none;fillColor=none;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;" vertex="1" parent="4owkeiqQ8d7qHkTqZ-NT-1">
+          <mxGeometry y="98" width="190" height="26" as="geometry" />
+        </mxCell>
+        <mxCell id="4owkeiqQ8d7qHkTqZ-NT-5" value="Extends" style="endArrow=block;endSize=16;endFill=0;html=1;entryX=1;entryY=1;entryDx=0;entryDy=0;exitX=0.5;exitY=0;exitDx=0;exitDy=0;" edge="1" parent="1" source="4owkeiqQ8d7qHkTqZ-NT-1" target="zvqzcq4YEUCvtmfPgBdI-33">
+          <mxGeometry y="10" width="160" relative="1" as="geometry">
+            <mxPoint x="514.7802197802198" y="490" as="sourcePoint" />
+            <mxPoint x="495" y="370" as="targetPoint" />
+            <mxPoint as="offset" />
+          </mxGeometry>
+        </mxCell>
+      </root>
+    </mxGraphModel>
+  </diagram>
+</mxfile>

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}];

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}]);
 

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

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);
+
+

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

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