adding plot_binned_strains_vs_radius

__init__.py

@@ -12,6 +12,7 @@ from plot_displacement_error import *
 from compute_displacements import *
 from compute_unwarped_images import *
 from image_expressions_py import *
+from plot_binned_strains_vs_radius import *
 from compute_displacement_error import *
 from compute_warped_images import *
 from compute_quadrature_degree import *

compute_strains.py

@@ -32,7 +32,8 @@ def compute_strains(
         CYL_or_PPS="PPS",
         write_strains=1,
         plot_strains=1,
-        write_strain_vs_radius=0,
+        write_strains_vs_radius=0,
+        write_binned_strains_vs_radius=0,
         verbose=1):
 
     if (mesh_w_local_basis_folder is not None) and (mesh_w_local_basis_basename is not None):
@@ -70,9 +71,13 @@ def compute_strains(
         strain_file = open(sol_folder+"/"+sol_basename+"-strains.dat", "w")
         strain_file.write("#t Err_avg Err_std Ecc_avg Ecc_std Ell_avg Ell_std Erc_avg Erc_std Erl_avg Erl_std Ecl_avg Ecl_std\n")
 
-    if (write_strain_vs_radius):
+    if (write_strains_vs_radius):
         strain_vs_radius_file = open(sol_folder+"/"+sol_basename+"-strains_vs_radius.dat", "w")
-        strain_vs_radius_file.write("#t r Err Ecc Ell Erc Erl Ecl\n")
+        strain_vs_radius_file.write("#t rr Err Ecc Ell Erc Erl Ecl\n")
+
+    if (write_binned_strains_vs_radius):
+        binned_strain_vs_radius_file = open(sol_folder+"/"+sol_basename+"-binned_strains_vs_radius.dat", "w")
+        binned_strain_vs_radius_file.write("#t rr Err Ecc Ell Erc Erl Ecl\n")
 
     if (ref_frame is not None):
         ref_mesh_filename = sol_folder+"/"+sol_basename+"_"+str(ref_frame).zfill(sol_zfill)+"."+sol_ext
@@ -118,7 +123,7 @@ def compute_strains(
             filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext,
             verbose=verbose)
 
-        if (write_strains) or (write_strain_vs_radius):
+        if (write_strains) or (write_strains_vs_radius) or (write_binned_strains_vs_radius):
             if (mesh_w_local_basis is not None):
                 assert (mesh.GetCellData().HasArray(strain_array_name+"_"+CYL_or_PPS))
                 farray_strain = mesh.GetCellData().GetArray(strain_array_name+"_"+CYL_or_PPS)
@@ -156,14 +161,35 @@ def compute_strains(
                     strain_file.write("".join([" " + str(strains_per_sector_avg[k_comp]) + " " + str(strains_per_sector_std[k_comp]) for k_comp in xrange(6)]))
             strain_file.write("\n")
 
-        if (write_strain_vs_radius):
-            assert (mesh_w_local_basis.GetCellData().HasArray("r"))
-            farray_r = mesh_w_local_basis.GetCellData().GetArray("r")
+        if (write_strains_vs_radius):
+            assert (mesh_w_local_basis.GetCellData().HasArray("rr"))
+            farray_rr = mesh_w_local_basis.GetCellData().GetArray("rr")
             for k_cell in xrange(n_cells):
-                strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, farray_r.GetTuple1(k_cell)]+list(farray_strain.GetTuple(k_cell))]) + "\n")
+                strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, farray_rr.GetTuple1(k_cell)]+list(farray_strain.GetTuple(k_cell))]) + "\n")
             strain_vs_radius_file.write("\n")
             strain_vs_radius_file.write("\n")
 
+        if (write_binned_strains_vs_radius):
+            assert (mesh_w_local_basis.GetCellData().HasArray("rr"))
+            farray_rr = mesh_w_local_basis.GetCellData().GetArray("rr")
+            n_r = 10
+            binned_strains = [[] for k_r in xrange(n_r)]
+            for k_cell in xrange(n_cells):
+                k_r = int(farray_rr.GetTuple1(k_cell)*n_r)
+                binned_strains[k_r].append(list(farray_strain.GetTuple(k_cell)))
+            #print binned_strains
+            binned_strains_avg = []
+            binned_strains_std = []
+            for k_r in xrange(n_r):
+                binned_strains_avg.append(numpy.mean(binned_strains[k_r], 0))
+                binned_strains_std.append(numpy.std (binned_strains[k_r], 0))
+            #print binned_strains_avg
+            #print binned_strains_std
+            for k_r in xrange(n_r):
+                binned_strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, (k_r+0.5)/n_r]+[val for k_comp in xrange(6) for val in [binned_strains_avg[k_r][k_comp], binned_strains_std[k_r][k_comp]]]]) + "\n")
+            binned_strain_vs_radius_file.write("\n")
+            binned_strain_vs_radius_file.write("\n")
+
     if (write_strains):
         strain_file.close()
 
@@ -174,5 +200,8 @@ def compute_strains(
                 suffix=None,
                 verbose=verbose)
 
-    if (write_strain_vs_radius):
+    if (write_strains_vs_radius):
         strain_vs_radius_file.close()
+
+    if (write_binned_strains_vs_radius):
+        binned_strain_vs_radius_file.close()

plot_binned_strains_vs_radius.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import math
+import numpy
+import os
+import sys
+
+########################################################################
+
+def plot_binned_strains_vs_radius(
+        working_folder,
+        working_basenames,
+        ref_folder=None,
+        ref_basename=None,
+        components="all",
+        suffix="",
+        verbose=1):
+
+    if (ref_folder is not None) and (ref_basename is not None):
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
+    else:
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()
+    n_frames = len(lines)-1
+
+    assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+    elif (components == "circ-long"):
+        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+    elif (components == "rad-circ"):
+        comp_names = ["radial", "circumferential", "radial-circumferential"]
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basenames[0]+"-binned_strains_vs_radius"
+    else:
+        plotfile_basename = "plot_binned_strains_vs_radius"+("-"+suffix)*(suffix!="")
+
+    plotfile = open(plotfile_basename+".plt", "w")
+    plotfile.write('''\
+set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
+
+set output "'''+plotfile_basename+'''.pdf"
+
+load "Set1.plt"
+set linestyle 1 pointtype 1
+set linestyle 2 pointtype 1
+set linestyle 3 pointtype 1
+set linestyle 4 pointtype 1
+set linestyle 5 pointtype 1
+set linestyle 6 pointtype 1
+set linestyle 7 pointtype 1
+set linestyle 8 pointtype 1
+set linestyle 9 pointtype 1
+
+set key box textcolor variable width +0
+
+set grid
+
+''')
+    for k_frame in xrange(n_frames):
+        plotfile.write('''\
+set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
+
+set xrange [0.:1.]
+set xtics("endocardium" 0, "epicardium" 1)
+
+''')
+        for comp_name in comp_names:
+            if   (comp_name == "radial"                      ): k_comp = 0
+            elif (comp_name == "circumferential"             ): k_comp = 1
+            elif (comp_name == "longitudinal"                ): k_comp = 2
+            elif (comp_name == "radial-circumferential"      ): k_comp = 3
+            elif (comp_name == "radial-longitudinal"         ): k_comp = 4
+            elif (comp_name == "circumferential-longitudinal"): k_comp = 5
+            plotfile.write('''\
+set ylabel "'''+comp_name+''' strain (%)"
+set yrange [-50:50]
+
+plot 0 linecolor rgb "black" notitle,\\
+''')
+            if (ref_folder is not None) and (ref_basename is not None):
+                plotfile.write('''\
+     "'''+ref_folder+'''/'''+ref_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+''') index '''+str(k_frame)+''' with lines linecolor "black" linewidth 3 notitle,\\
+     "'''+ref_folder+'''/'''+ref_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+'''):(100*$'''+str(3+2*k_comp+1)+''') index '''+str(k_frame)+''' with errorbars linecolor "black" linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
+''')+(len(working_basenames)==0)*('''
+
+'''))
+            for k_basename in range(len(working_basenames)):
+                working_basename = working_basenames[k_basename]
+                plotfile.write('''\
+     "'''+working_folder+'''/'''+working_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+''') index '''+str(k_frame)+''' with lines linestyle '''+str(k_basename+1)+''' linewidth 3 title "'''+working_basename+'''",\\
+     "'''+working_folder+'''/'''+working_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+'''):(100*$'''+str(3+2*k_comp+1)+''') index '''+str(k_frame)+''' with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
+''')+(k_basename==len(working_basenames)-1)*('''
+
+'''))
+
+    plotfile.close()
+
+    os.system("gnuplot "+plotfile_basename+".plt")

plot_strains.py

@@ -22,11 +22,10 @@ def plot_strains(
         suffix="",
         verbose=1):
 
-
-    if (working_scalings is None):
-        working_scalings = [1.] * len(working_basenames)
-    else:
+    if (working_scalings is not None):
         assert (len(working_scalings) == len(working_basenames))
+    else:
+        working_scalings = [1.] * len(working_basenames)
 
     if (ref_folder is not None) and (ref_basename is not None):
         lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
@@ -38,12 +37,19 @@ def plot_strains(
     #print "n_sectors = " + str(n_sectors)
 
     assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+    elif (components == "circ-long"):
+        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+    elif (components == "rad-circ"):
+        comp_names = ["radial", "circumferential", "radial-circumferential"]
 
     if (suffix is None):
         plotfile_basename = working_folder+"/"+working_basenames[0]+"-strains"
     else:
         plotfile_basename = "plot_strains"+("-"+suffix)*(suffix!="")
     plotfile = open(plotfile_basename+".plt", "w")
+
     plotfile.write('''\
 set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
 
@@ -70,15 +76,10 @@ set grid
         plotfile.write('''\
 set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
 
+set xlabel "frame ()"
+set xrange [0:'''+str(n_frames)+''']
+
 ''')
-        if (components == "all"):
-            comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
-        elif (components == "circ-long"):
-            comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
-        elif (components == "rad-circ"):
-            comp_names = ["radial", "circumferential", "radial-circumferential"]
-        else:
-            assert (0), "components must be \"all\", \"circ-long\" or \"rad-circ\". Aborting."
         for comp_name in comp_names:
             if   (comp_name == "radial"                      ): k_comp = 0
             elif (comp_name == "circumferential"             ): k_comp = 1
@@ -87,11 +88,18 @@ set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in
             elif (comp_name == "radial-longitudinal"         ): k_comp = 4
             elif (comp_name == "circumferential-longitudinal"): k_comp = 5
             plotfile.write('''\
-set xrange [0:'''+str(n_frames)+''']
-
 set ylabel "'''+comp_name+''' strain (%)"
+''')
+            if ("-" in comp_name):
+                plotfile.write('''\
+set yrange [-10:10]
+''')
+            else:
+                plotfile.write('''\
 set yrange [-30:30]
 
+''')
+            plotfile.write('''\
 plot 0 linecolor rgb "black" notitle,\\
 ''')
             if (ref_folder is not None) and (ref_basename is not None):
@@ -105,8 +113,8 @@ plot 0 linecolor rgb "black" notitle,\\
                 working_basename = working_basenames[k_basename]
                 working_scaling = working_scalings[k_basename]
                 plotfile.write('''\
-     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./100):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\
-     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./100):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
+     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\
+     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
 ''')+(k_basename==len(working_basenames)-1)*('''
 
 '''))

plot_strains_vs_radius.py

@@ -18,20 +18,36 @@ import sys
 def plot_strains_vs_radius(
         working_folder,
         working_basenames,
-        index=10,
         ref_folder=None,
         ref_basename=None,
         components="all",
         suffix="",
         verbose=1):
 
+    if (ref_folder is not None) and (ref_basename is not None):
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
+    else:
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()
+    n_frames = len(lines)-1
+
     assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+    elif (components == "circ-long"):
+        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+    elif (components == "rad-circ"):
+        comp_names = ["radial", "circumferential", "radial-circumferential"]
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basenames[0]+"-strains_vs_radius"
+    else:
+        plotfile_basename = "plot_strains_vs_radius"+("-"+suffix)*(suffix!="")
 
-    plotfile = open("plot_strains_vs_radius"+("-"+suffix)*(suffix!="")+".plt", "w")
+    plotfile = open(plotfile_basename+".plt", "w")
     plotfile.write('''\
 set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
 
-set output "plot_strains_vs_radius'''+('''-'''+suffix)*(suffix!="")+'''.pdf"
+set output "'''+plotfile_basename+'''.pdf"
 
 load "Set1.plt"
 set linestyle 1 pointtype 1
@@ -48,48 +64,42 @@ set key box textcolor variable width +0
 
 set grid
 
-set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
-
 ''')
-    if (components == "all"):
-        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
-    elif (components == "circ-long"):
-        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
-    elif (components == "rad-circ"):
-        comp_names = ["radial", "circumferential", "radial-circumferential"]
-    else:
-        assert (0), "components must be \"all\", \"circ-long\" or \"rad-circ\". Aborting."
-    for comp_name in comp_names:
-        if   (comp_name == "radial"                      ): k_comp = 0
-        elif (comp_name == "circumferential"             ): k_comp = 1
-        elif (comp_name == "longitudinal"                ): k_comp = 2
-        elif (comp_name == "radial-circumferential"      ): k_comp = 3
-        elif (comp_name == "radial-longitudinal"         ): k_comp = 4
-        elif (comp_name == "circumferential-longitudinal"): k_comp = 5
+    for k_frame in xrange(n_frames):
         plotfile.write('''\
+set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
 
-set xlabel "r (mm)"
-set xrange [0.2:0.4]
+set xrange [0.:1.]
+set xtics("endocardium" 0, "epicardium" 1)
 
+''')
+        for comp_name in comp_names:
+            if   (comp_name == "radial"                      ): k_comp = 0
+            elif (comp_name == "circumferential"             ): k_comp = 1
+            elif (comp_name == "longitudinal"                ): k_comp = 2
+            elif (comp_name == "radial-circumferential"      ): k_comp = 3
+            elif (comp_name == "radial-longitudinal"         ): k_comp = 4
+            elif (comp_name == "circumferential-longitudinal"): k_comp = 5
+            plotfile.write('''\
 set ylabel "'''+comp_name+''' strain (%)"
-set yrange [-40:40]
+set yrange [-50:50]
 
 plot 0 linecolor rgb "black" notitle,\\
 ''')
-        if (ref_folder is not None) and (ref_basename is not None):
-            plotfile.write('''\
-     "'''+ref_folder+'''/'''+ref_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(index)+''' with points linecolor "black" pointtype 1 pointsize 1 linewidth 3 notitle'''+(len(working_basenames)>0)*(''',\\
+            if (ref_folder is not None) and (ref_basename is not None):
+                plotfile.write('''\
+     "'''+ref_folder+'''/'''+ref_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(k_frame)+''' with points linecolor "black" pointtype 1 pointsize 1 linewidth 3 notitle'''+(len(working_basenames)>0)*(''',\\
 ''')+(len(working_basenames)==0)*('''
 
 '''))
-        for k_basename in range(len(working_basenames)):
-            working_basename = working_basenames[k_basename]
-            plotfile.write('''\
-     "'''+working_folder+'''/'''+working_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(index)+''' with points linestyle '''+str(k_basename+1)+''' pointtype 1 pointsize 1 linewidth 3 title "'''+working_basename+'''"'''+(k_basename<len(working_basenames)-1)*(''',\\
+            for k_basename in range(len(working_basenames)):
+                working_basename = working_basenames[k_basename]
+                plotfile.write('''\
+     "'''+working_folder+'''/'''+working_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(k_frame)+''' with points linestyle '''+str(k_basename+1)+''' pointtype 1 pointsize 1 linewidth 3 title "'''+working_basename+'''"'''+(k_basename<len(working_basenames)-1)*(''',\\
 ''')+(k_basename==len(working_basenames)-1)*('''
 
 '''))
 
     plotfile.close()
 
-    os.system("gnuplot plot_strains_vs_radius"+("-"+suffix)*(suffix!="")+".plt")
+    os.system("gnuplot "+plotfile_basename+".plt")