Calculate PCs

main.js

@@ -5,6 +5,65 @@ async function main() {
       nv1.setColormap(nv1.volumes[2].id, 'red')
       lesionSlider.oninput()
   }
+  function norm0to1(val, mn, mx) {
+    // return clamped/normalized 0..1, linearly interpolated min..max
+    let range = mx - mn;
+    let ret = (val - mn)/range;
+    return Math.min(Math.max(ret,0),1);
+  }
+  function neglect_predict() {
+    if (nv1.volumes.length !== 3) {
+      window.alert('Please reload the page and open lesions with the "Open Lesion Map" button')
+      return
+    }
+    let mask = nv1.volumes[1].img
+    let lesion = nv1.volumes[2].img
+    let nvox = lesion.length
+    if (nvox !== mask.length) {
+      window.alert('Lesion must precisely match mask')
+      return
+    }
+    let lesionVol = 0
+    let nMask = 0
+    const nVoxPCA = 5
+    const nMaskX = Math.floor(pca_coeff.length / nVoxPCA) //expected size
+    let map = new Float64Array(nvox)
+    let v = 0
+    for (let i = 0; i < nvox; i++) {
+      if ((lesion[i] > 0) && (mask[i] > 0))
+        lesionVol ++
+      if (mask[i] > 1)
+        map[nMask++] = lesion[i]
+    }
+    if (nMaskX !== nMask) {
+      window.alert('PCA and mask have different number of elements')
+      return
+    }
+    let PC = new Float64Array(nVoxPCA)
+    //PC = (map'-pca_val.mu)*pca_val.coeff
+    for (let i = 0; i < nMask; i++) {
+      map[i] -= pca_mu[i]
+    }
+    for (let p = 0; p < nVoxPCA; p++) {
+      const offset = p * nMask
+      for (let i = 0; i < nMask; i++) {
+        PC[p] += map[i]*pca_coeff[i+offset]
+      }
+    }
+    //normalize values to range 0..1
+    // -> PCs: min = -51.9073; max = 110.0535
+    // -> CoC: min = -0.024243014; max = 0.951938077
+    // -> ROI_vol: min = 0; max = 21.625
+    for (let p = 0; p < nVoxPCA; p++) {
+        PC[p] = norm0to1(PC[p], -51.9073, 110.0535)
+    }
+    let acuteCoC = parseFloat(cocNumber.value)
+    acuteCoC = norm0to1(acuteCoC, -0.024243014, 0.951938077)
+    let ROI_volML = lesionVol / 1000
+    ROI_volML = norm0to1(ROI_volML, 0, 21.625);
+    const input_vector = [PC[0], PC[1], PC[2], PC[3], PC[4], acuteCoC, ROI_volML];
+    console.log(input_vector)
+  }
   openBtn.onclick = async function () {
     let input = document.createElement('input')
     input.style.display = 'none'
@@ -39,18 +98,21 @@ async function main() {
   }
   const nv1 = new Niivue(defaults)
   nv1.attachToCanvas(gl1)
+  const pca_coeff = (await nv1.loadFromUrl('./pca_values_coeff.nii.gz')).img
+  const pca_mu = (await nv1.loadFromUrl('./pca_values_mu.nii.gz')).img
   nv1.opts.dragMode = nv1.dragModes.pan
   nv1.opts.multiplanarForceRender = true
   nv1.opts.yoke3Dto2DZoom = true
   nv1.opts.crosshairGap = 11
   nv1.setInterpolation(true)
   await nv1.loadVolumes([
     { url: './betsct1_unsmooth.nii.gz' },
-    {url: './mask_vox.nii.gz', colormap: "blue"},
+    {url: './mask.nii.gz', colormap: "blue"},
     {url: './M2095_lesion.nii.gz', colormap: "red"},
   ])
   maskSlider.oninput()
   lesionSlider.oninput()
+  neglect_predict()
 }
 
 main()

public/neglect_predict.m

@@ -0,0 +1,100 @@
+function neglect_predict(fnm, acuteCoC)
+%Predict chronic recovery based on lesion map and acute center of
+%cancellation score.
+%Examples
+% neglect_predict; %use GUI;
+% neglect_predict('M2095_lesion.nii.gz', 0.65)
+    if nargin < 1
+        [p_file, p_path] = uigetfile('*.nii.gz;*.nii', 'Select lesion map');
+        if p_file==0
+            return
+        end
+        fnm = fullfile(p_path, p_file);
+    end
+    if nargin < 2
+        prompt = {'Enter acute CoC score (-1..1):'};
+        dlgtitle = 'Input';
+        fieldsize = [1];
+        definput = {'0.65'};
+        answer = inputdlg(prompt,dlgtitle,fieldsize,definput);
+        if length(answer) == 0
+            return
+        end
+        acuteCoC = str2double(answer{1});
+    end
+    lesion = niftiread(fnm);
+    mpath = fileparts(mfilename("fullpath"));
+    fnmMsk = fullfile(mpath, 'mask.nii.gz');
+    if ~exist(fnmMsk,'file')
+        error('Unable to find %s', fnmMsk)
+    end
+    maskVox012 = uint8(niftiread(fnmMsk));
+    maskROI = uint8(maskVox012 > 0);
+    maskVox = uint8(maskVox012 > 1);
+    ROI_volVox = nnz(maskROI & lesion);
+    ROI_volML = ROI_volVox / 1000; %convert voxels to ML
+    fprintf("%d lesioned voxels in ROI mask: %g ml\n", ROI_volVox, ROI_volML);
+    %PCA
+    fnmPCA = fullfile(mpath, 'pca_values_5x21220.mat');
+    if ~exist(fnmPCA,'file')
+        error('Unable to find %s', fnmPCA)
+    end
+    pca_val = load(fnmPCA).pca_val;
+    if size(pca_val.mu,2) ~= nnz(maskVox)
+        error('maskVox size does not match pca_val')
+    end
+    map_org = double(lesion(:));
+    mask_log = maskVox ~=0;
+    map = map_org(mask_log);
+    PC = (map'-pca_val.mu)*pca_val.coeff;
+    if numel(PC) ~= 5
+        error('scores_new should have 5 values')
+    end
+    %normalize values to range 0..1
+    %       -> PCs: min = -51.9073; max = 110.0535
+    %       -> CoC: min = -0.024243014; max = 0.951938077
+    %       -> ROI_vol: min = 0; max = 21.625
+    for i = 1:5
+        PC(i) = norm0to1(PC(i), -51.9073, 110.0535);
+    end
+   acuteCoC = norm0to1(acuteCoC, -0.024243014, 0.951938077);
+    ROI_volML = norm0to1(ROI_volML, 0, 21.625);
+    % input_vector = [PC1, PC2, PC3, PC4, PC5, CoC, ROI_vol];
+    % Exemplary patient data (with outcome = 0.8155):
+    input_vector = [PC(1), PC(2), PC(3), PC(4), PC(5), acuteCoC, ROI_volML];
+    fnmModel = fullfile(mpath, 'models_5x10_diff.mat');
+    if ~exist(fnmPCA,'file')
+        error('Unable to find %s', fnmPCA)
+    end
+    models = load(fnmModel).models;
+    % Blank prediction-array
+    predictions_mdls = zeros(numel(models),1);
+
+    % For each single model (5-fold nested cross validation x 10 model repetitions)
+    for i = 1:numel(models)
+        % Extract information from the i-th model
+        support_vectors_i = models{i}.SVs;
+        coefficients_i = models{i}.sv_coef;
+        bias_i = -models{i}.rho;
+        % Define RBF kernel function
+        gamma_i = models{i}.Parameters(4);
+        rbf_kernel = @(x1, x2) exp(-gamma_i * sum((x1 - x2).^2));
+        % Calculate the kernel values
+        kernel_values_i = arrayfun(@(j) rbf_kernel(input_vector, support_vectors_i(j, :)), 1:size(support_vectors_i, 1));
+        % Calculate the prediction using the regression function
+        predictions_mdls(i) = sum(coefficients_i' .* kernel_values_i) + bias_i;
+        % Feature weights and bias term
+        % w = coefficients_i' * support_vectors_i; 
+        % b = bias_i;
+    end
+    % Calculate mean prediction
+    prediction_mean = mean(predictions_mdls);
+    disp(['Mean prediction: ' num2str(prediction_mean)]);
+end
+function ret = norm0to1(val, mn, mx)
+    %return normalized 0..1, linearly interpolated min..max
+    range = mx - mn;
+    ret = (val - mn)/range;
+    ret = min(max(ret,0),1);
+end
+