A new version of the propagation velocity script have been implemented

scripts/calc_propagation_velocity_full_grid/CMakeLists.txt

@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 2.8)
+
+PROJECT(CalcPropagationVelocity)
+
+find_package(VTK REQUIRED)
+include(${VTK_USE_FILE})
+
+SET( CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin )
+
+add_executable(CalcPropagationVelocity main.cpp )
+
+target_link_libraries(CalcPropagationVelocity ${VTK_LIBRARIES})

scripts/calc_propagation_velocity_full_grid/calc_activation_time_full_grid.py

@@ -0,0 +1,181 @@
+# Author: Lucas Berg
+#
+# Program that calculates the APD of every cell in the grid by passing a certain percentage. 
+#   e.g: APD_90 --> percentage = 90
+#        APD_80 --> percentage = 80
+#        APD_70 --> percentage = 70
+#
+# This program also works with multiples AP's !
+
+import sys
+import subprocess
+import time
+import numpy as np
+
+def forwarddiff(y, h):
+    n = len(y)
+    res = []
+
+    for i in range(1, n):
+        res.append((y[i] - y[i-1]) / h)
+
+    return res
+
+
+def slope_start(data, start=0, epsilon=0.0001, h=1.0):
+
+    d = data[start:]
+    n = len(d)
+
+    for i in range(1,n):
+        if abs(d[i] - d[i-1]/h) > epsilon:
+            return i+start
+
+
+def slope_end(data, start=0, epsilon=0.0001, h=1.0):
+
+    d = data[start:]
+    n = len(d)
+
+    for i in range(1,n):
+        if abs(d[i] - d[i-1]/h) < epsilon:
+            return i+start
+
+
+def max_index(data, start, end):
+
+    d = data[start:(start+end)]
+
+    max_v = max(d)
+
+    max_index = d.index(max_v) + start
+
+    return max_index
+
+
+def index_activation(data, start=0):
+
+    d = data[start:]
+
+    for i, v in enumerate(d):
+        if d[i + start] < 0.0 < d[i + start + 1]:
+            return i+start
+
+def calc_ap_limits (vms, start, end):
+
+    v = vms[start:end].tolist()
+
+    maxV = max(v)
+    minV = min(v)
+
+    index_maxV = v.index(maxV) + start
+
+    return index_maxV, maxV, minV
+
+def calc_reference_potential (minV, maxV, percentage):
+
+    refV = minV + (maxV - minV)*(1.0 - percentage)
+
+    return refV        
+
+def calc_time_reference (vms, start, end, h, refV):
+
+    v = vms[start:(start+end)]
+
+    for i in range(len(v)):
+        if (v[i] < refV):
+            return (i+start), (i+start)*h
+    return -1, -1
+
+def calc_max_min_ref_potential (data,apd_p,start,end):
+
+    d = data[start:end]
+
+    max_v = max(d)
+    min_v = min(d)
+    ref_v = min_v + (max_v - min_v)*(1.0 - apd_p)
+
+    return max_v, min_v, ref_v
+
+def calc_max_derivative (vms,start,end,h):
+    v = vms[start:end]
+    n = len(v)
+    dvdt = range(0,n)
+
+    for i in range(1,n):
+        dvdt[i] = abs(v[i] - v[i-1]/h)
+
+    max_dvdt = max(dvdt)
+
+    max_index = dvdt.index(max_dvdt) + start
+
+    return max_index, max_index*h
+
+def calc_activation_time (vms, h, start, end):
+
+    index_peak, t_peak = calc_max_derivative(vms,start,end,h)
+
+    return (t_peak)
+
+def main ():
+    if ( len(sys.argv) != 5 ):
+        print("-------------------------------------------------------------------------------------------------------------")
+        print("Usage:> python %s <output_dir> <ms_each_step> <total_num_cells> <print_rate>" % sys.argv[0])
+        print("-------------------------------------------------------------------------------------------------------------")
+        print("<output_dir> = Output directory of the simulation")
+        print("<ms_each_step> = Number of milliseconds from each timestep")
+        print("     this value can be calculated as:")
+        print("         num_files = (simulation_time) / (dt * print_rate)")
+        print("         ms_each_step = (simulation_time) / (num_files)")
+        print("     Where the values <simulation_time>, <dt> and <print_rate> are all")
+        print("     given in the configuration file of the simulation.")
+        print("<total_num_cells> = Total number of cells to calculate the APD")
+        print("<print_rate> = The print rate of the simulation")
+        print("-------------------------------------------------------------------------------------------------------------")
+        print("Example:> python calc_activation_time_full_grid.py ../../outputs/plain_100_100_100_fhn 2 10000 100")
+        print("          python calc_activation_time_full_grid.py ../../outputs/plain_100_100_100_tentusscher 2 10000 100")
+        print("-------------------------------------------------------------------------------------------------------------")
+        return 1
+
+    # Get user inputs
+    output_dir = sys.argv[1]
+    ms_each_step = float(sys.argv[2])
+    total_num_cells = int(sys.argv[3])
+    print_rate = int(sys.argv[4])
+
+    # Get the transmembrane potential for all timesteps and every single cell in the grid
+    print("[!] Calling 'getAps.sh' script ...")
+    cmd = "./getAps.sh %s V 6 %d %d" % (output_dir,total_num_cells,print_rate)
+    rc = subprocess.call( cmd, shell=True )
+
+    # Open the generated file from the previous script as a Numpy array and get its dimensions
+    timesteps_file = open("timesteps.txt")
+    ap_data = np.genfromtxt(timesteps_file)
+    num_cells = ap_data.shape[0]
+    num_timesteps = ap_data.shape[1]
+
+    # Open the output file
+    out_file = open("inputs/cells-at.txt","w")
+
+    print("[!] Calculating Activation Times ...")
+
+    begin = time.time()
+    # Iterate over each cell
+    for i in range(num_cells):
+        # Get the transmembrane potential from the current cell
+        vms = ap_data[i]
+
+        # Calculate the activation time of that cell
+        at = calc_activation_time(vms,ms_each_step,0,500)
+
+        # Write the value to an output file
+        out_file.write("%g\n" % (at))
+    end = time.time()
+
+    print("[!] Elapsed time = %g seconds" % (end - begin))
+
+    out_file.close()
+    print("[+] Output file saved in 'inputs/cell-at.txt'")
+
+if __name__ == "__main__":
+    main()

scripts/calc_propagation_velocity_full_grid/clean_files.sh

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm timesteps.txt
+rm aps/*.txt

scripts/calc_propagation_velocity_full_grid/clean_project.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+BUILD_DIR="build"
+
+if [[ "$#" -eq 1 ]]; then
+    BUILD_DIR=$1
+fi
+
+if [[ ! -d "${BUILD_DIR}" ]]; then
+  echo "Directory ${BUILD_DIR} does not exist"
+  exit
+fi
+
+rm -fr ${BUILD_DIR}/*
+rm -fr bin/*

scripts/calc_propagation_velocity_full_grid/getAps.sh

@@ -0,0 +1,46 @@
+#!/usr/bin/env bash
+AP_DIR=$1
+AP_PREFIX=$2
+AP_LINE=$3
+TOTAL_NUMBER_CELLS=$4
+PRINT_RATE=$5
+
+if [[ "$#" -ne 5 ]]; then
+    echo "---------------------------------------------------------------------------------"
+    echo "Usage:> $0 <AP_DIR> <AP_PREFIX> <AP_LINE> <TOTAL_NUMBER_CELLS> <PRINT_RATE>"
+    echo "---------------------------------------------------------------------------------"
+    echo "<AP_DIR> = Directory where the results of the simulation are stored"
+    echo "<AP_PREFIX> = Prefix of the results files"
+    echo "<AP_LINE> = Line where the action potential information is stored"
+    echo "<TOTAL_NUMBER_CELLS> = Total number of grid cells"
+    echo "<PRINT_RATE> = The print rate of the simulation"
+    echo "---------------------------------------------------------------------------------"
+    echo "!!! This script only works if the output file is NOT saved in   !!!"
+    echo "!!!                     binary format                           !!!"
+    echo "---------------------------------------------------------------------------------"
+    echo "<AP_LINE> guide:"
+    echo "  VTP/VTU --> Line 6"
+    echo "  VTK     --> Last line"
+    echo "---------------------------------------------------------------------------------"
+    exit 1
+fi
+
+# Get the transmembrane potential from each of the .vtu files and dump it into a .txt
+TIMESTEP=0
+for i in `ls -1v ${AP_DIR}/${AP_PREFIX}*`; do
+    AP_OUT="aps/timestep-$TIMESTEP.txt"
+    sed -n "${AP_LINE}p" $i | awk -v TOTAL_NUMBER_CELLS="$TOTAL_NUMBER_CELLS" -F ' ' '{ for (i=0; i < TOTAL_NUMBER_CELLS; i++) { printf "%g\n", $i } }' > ${AP_OUT}
+    let "TIMESTEP=TIMESTEP+$PRINT_RATE"
+done
+
+# Adjust the total number of timesteps
+let "TOTAL_TIMESTEPS=TIMESTEP-$PRINT_RATE"
+
+# Concatenate the filename from each timestep
+CMD=""
+for i in $(seq 0 $PRINT_RATE $TOTAL_TIMESTEPS); do
+    CMD="$CMD aps/timestep-$i.txt"
+done
+
+# Use the 'paste' command to append each timestep into a column on the output file
+paste $CMD > timesteps.txt