--- a
+++ b/code/meshfitting.py
@@ -0,0 +1,492 @@
+import os
+from multiprocessing import Pool
+from functools import partial
+
+
+def meshGeneration(subject_dir, template_dir, param_dir, tmps_dir, vtks_dir, dofs_dir): 
+   
+    os.system('rm {0}/*.txt'.format(subject_dir))  
+    
+    os.system('rm {0}/*'.format(vtks_dir))    
+    
+    os.system('rm {0}/*'.format(tmps_dir))    
+
+    for fr in ['ED', 'ES']:
+        
+        ###########################################################################    
+        # extract meshes of lvendo, lvepi, lvmyo, rv and rveip at ED and ES, respectively
+        os.system('binarize '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_RV_{2}.nii.gz '
+                  '4 4 255 0'
+                  .format(subject_dir, tmps_dir, fr))
+           
+        os.system('mcubes '
+                  '{0}/vtk_RV_{2}.nii.gz '
+                  '{1}/RV_{2}.vtk '
+                  '120 -blur 2'
+                  .format(tmps_dir, vtks_dir, fr))
+        
+        os.system('binarize '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_RVepi_{2}.nii.gz '
+                  '3 4 255 0'
+                  .format(subject_dir, tmps_dir, fr))
+        
+        os.system('mcubes '
+                  '{0}/vtk_RVepi_{2}.nii.gz '
+                  '{1}/RVepi_{2}.vtk '
+                  '120 -blur 2'
+                  .format(tmps_dir, vtks_dir, fr))
+        
+        os.system('padding '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_LV_{2}.nii.gz '
+                  '4 0'
+                  .format(subject_dir, tmps_dir, fr))
+        
+        os.system('padding '
+                  '{1}/vtk_LV_{2}.nii.gz '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_LV_{2}.nii.gz '
+                  '3 0'
+                  .format(subject_dir, tmps_dir, fr))
+        
+        os.system('binarize '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_LVendo_{2}.nii.gz '
+                  '1 1 255 0'
+                  .format(subject_dir, tmps_dir, fr))
+        
+        os.system('mcubes '
+                  '{0}/vtk_LVendo_{2}.nii.gz '
+                  '{1}/LVendo_{2}.vtk '
+                  '120 -blur 2'
+                  .format(tmps_dir, vtks_dir, fr))
+        
+        os.system('binarize '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_LVepi_{2}.nii.gz '
+                  '1 2 255 0'
+                  .format(subject_dir, tmps_dir, fr))
+        
+        os.system('mcubes '
+                  '{0}/vtk_LVepi_{2}.nii.gz '
+                  '{1}/LVepi_{2}.vtk '
+                  '120 -blur 2'
+                  .format(tmps_dir, vtks_dir, fr))
+        
+        os.system('binarize '
+                  '{0}/PHsegmentation_{2}.gipl '
+                  '{1}/vtk_LVmyo_{2}.nii.gz '
+                  '2 2 255 0'
+                  .format(subject_dir, tmps_dir, fr))
+        
+        os.system('mcubes '
+                  '{0}/vtk_LVmyo_{2}.nii.gz '
+                  '{1}/LVmyo_{2}.vtk '
+                  '120 -blur 2'
+                  .format(tmps_dir, vtks_dir, fr))
+        
+#        os.system('binarize '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_RV_{2}.nii.gz '
+#                  '4 4 255 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#           
+#        os.system('mcubes '
+#                  '{0}/vtk_RV_{2}.nii.gz '
+#                  '{1}/RV_{2}.vtk '
+#                  '120 -blur 2'
+#                  .format(tmps_dir, vtks_dir, fr))
+#        
+#        os.system('binarize '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_RVepi_{2}.nii.gz '
+#                  '3 4 255 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#        
+#        os.system('mcubes '
+#                  '{0}/vtk_RVepi_{2}.nii.gz '
+#                  '{1}/RVepi_{2}.vtk '
+#                  '120 -blur 2'
+#                  .format(tmps_dir, vtks_dir, fr))
+#        
+#        os.system('padding '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_LV_{2}.nii.gz '
+#                  '4 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#        
+#        os.system('padding '
+#                  '{1}/vtk_LV_{2}.nii.gz '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_LV_{2}.nii.gz '
+#                  '3 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#        
+#        os.system('binarize '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_LVendo_{2}.nii.gz '
+#                  '1 1 255 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#        
+#        os.system('mcubes '
+#                  '{0}/vtk_LVendo_{2}.nii.gz '
+#                  '{1}/LVendo_{2}.vtk '
+#                  '120 -blur 2'
+#                  .format(tmps_dir, vtks_dir, fr))
+#        
+#        os.system('binarize '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_LVepi_{2}.nii.gz '
+#                  '1 2 255 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#        
+#        os.system('mcubes '
+#                  '{0}/vtk_LVepi_{2}.nii.gz '
+#                  '{1}/LVepi_{2}.vtk '
+#                  '120 -blur 2'
+#                  .format(tmps_dir, vtks_dir, fr))
+#        
+#        os.system('binarize '
+#                  '{0}/seg_lvsa_SR_{2}.nii.gz '
+#                  '{1}/vtk_LVmyo_{2}.nii.gz '
+#                  '2 2 255 0'
+#                  .format(subject_dir, tmps_dir, fr))
+#        
+#        os.system('mcubes '
+#                  '{0}/vtk_LVmyo_{2}.nii.gz '
+#                  '{1}/LVmyo_{2}.vtk '
+#                  '120 -blur 2'
+#                  .format(tmps_dir, vtks_dir, fr))
+    
+    ###############################################################################
+    #use landmark to initialise the registration
+#    os.system('prreg '
+#              '{0}/landmarks.vtk '
+#              '{1}/landmarks.vtk '
+#              '-dofout {2}/landmarks.dof.gz'
+#              .format(subject_dir, template_dir, dofs_dir))
+    
+    ###############################################################################
+    for fr in ['ED', 'ES']:
+    
+#        os.system('msrreg 3 '
+#                  '{0}/RV_{4}.vtk '
+#                  '{0}/LVendo_{4}.vtk '
+#                  '{0}/LVepi_{4}.vtk '
+#                  '{1}/RV_{4}.vtk '
+#                  '{1}/LVendo_{4}.vtk '
+#                  '{1}/LVepi_{4}.vtk '
+#                  '-dofin {2}/landmarks.dof.gz '
+#                  '-dofout {3}/{4}.dof.gz '
+#                  '-symmetric'
+#                  .format(vtks_dir, template_dir, dofs_dir, tmps_dir, fr))
+        
+        os.system('msrreg 3 '
+                  '{0}/RV_{3}.vtk '
+                  '{0}/LVendo_{3}.vtk '
+                  '{0}/LVepi_{3}.vtk '
+                  '{1}/RV_{3}.vtk '
+                  '{1}/LVendo_{3}.vtk '
+                  '{1}/LVepi_{3}.vtk '
+                  '-dofout {2}/{3}.dof.gz '
+                  '-symmetric'
+                  .format(vtks_dir, template_dir, tmps_dir, fr))
+        
+        os.system('msrreg 2 '
+                  '{0}/LVendo_{3}.vtk '
+                  '{0}/LVepi_{3}.vtk '
+                  '{1}/LVendo_{3}.vtk '
+                  '{1}/LVepi_{3}.vtk '
+                  '-dofin {2}/{3}.dof.gz '
+                  '-dofout {2}/lv_{3}_rreg.dof.gz '
+                  '-symmetric'
+                  .format(vtks_dir, template_dir, tmps_dir, fr))
+        
+        os.system('srreg '
+                  '{0}/RV_{3}.vtk '
+                  '{1}/RV_{3}.vtk '
+                  '-dofin {2}/{3}.dof.gz '
+                  '-dofout {2}/rv_{3}_rreg.dof.gz '
+                  '-symmetric'
+                  .format(vtks_dir, template_dir, tmps_dir, fr))
+        
+        ###########################################################################     
+        os.system('ptransformation '
+                  '{0}/RV_{2}.vtk '
+                  '{0}/N_RV_{2}.vtk '
+                  '-dofin {1}/rv_{2}_rreg.dof.gz'
+                  .format(vtks_dir, tmps_dir, fr))
+        
+        os.system('ptransformation '
+                  '{0}/RVepi_{2}.vtk '
+                  '{0}/N_RVepi_{2}.vtk '
+                  '-dofin {1}/rv_{2}_rreg.dof.gz'
+                  .format(vtks_dir, tmps_dir, fr))
+        
+        os.system('ptransformation '
+                  '{0}/LVendo_{2}.vtk '
+                  '{0}/N_LVendo_{2}.vtk '
+                  '-dofin {1}/lv_{2}_rreg.dof.gz'
+                  .format(vtks_dir, tmps_dir, fr))
+        
+        os.system('ptransformation '
+                  '{0}/LVepi_{2}.vtk '
+                  '{0}/N_LVepi_{2}.vtk '
+                  '-dofin {1}/lv_{2}_rreg.dof.gz'
+                  .format(vtks_dir, tmps_dir, fr))
+        
+        os.system('ptransformation '
+                  '{0}/LVmyo_{2}.vtk '
+                  '{0}/N_LVmyo_{2}.vtk '
+                  '-dofin {1}/lv_{2}_rreg.dof.gz'
+                  .format(vtks_dir, tmps_dir, fr))
+        
+        os.system('transformation '
+                  '{0}/vtk_RV_{1}.nii.gz '
+                  '{0}/N_vtk_RV_{1}.nii.gz '
+                  '-dofin {0}/lv_{1}_rreg.dof.gz '
+                  '-invert'
+                  .format(tmps_dir, fr))
+        
+        os.system('transformation '
+                  '{0}/vtk_LV_{1}.nii.gz '
+                  '{0}/N_vtk_LV_{1}.nii.gz '
+                  '-dofin {0}/lv_{1}_rreg.dof.gz '
+                  '-invert'
+                  .format(tmps_dir, fr))
+        
+        ###########################################################################
+        #affine
+        os.system('areg '
+                  '{0}/vtk_RV_{3}.nii.gz '
+                  '{1}/N_vtk_RV_{3}.nii.gz '
+                  '-dofout {1}/rv_{3}_areg.dof.gz '
+                  '-parin {2}/segareg.txt'
+                  .format(template_dir, tmps_dir, param_dir, fr))
+        
+        os.system('areg '
+                  '{0}/vtk_LV_{3}.nii.gz '
+                  '{1}/N_vtk_LV_{3}.nii.gz '
+                  '-dofout {1}/lv_{3}_areg.dof.gz '
+                  '-parin {2}/segareg.txt'
+                  .format(template_dir, tmps_dir, param_dir, fr))
+        
+        #non-rigid
+        os.system('nreg '
+                  '{0}/vtk_RV_{3}.nii.gz '
+                  '{1}/N_vtk_RV_{3}.nii.gz '
+                  '-dofin {1}/rv_{3}_areg.dof.gz '
+                  '-dofout {1}/rv_{3}_nreg.dof.gz '
+                  '-parin {2}/segreg.txt'
+                  .format(template_dir, tmps_dir, param_dir, fr))
+        
+        os.system('snreg '
+                  '{0}/RV_{3}.vtk '
+                  '{1}/N_RV_{3}.vtk '
+                  '-dofin {2}/rv_{3}_nreg.dof.gz '
+                  '-dofout {2}/rv{3}ds8.dof.gz '
+                  '-ds 8 -symmetric'
+                  .format(template_dir, vtks_dir, tmps_dir, fr))
+        
+        os.system('nreg '
+                  '{0}/vtk_LV_{3}.nii.gz '
+                  '{1}/N_vtk_LV_{3}.nii.gz '
+                  '-dofin {1}/lv_{3}_areg.dof.gz '
+                  '-dofout {1}/lv_{3}_nreg.dof.gz '
+                  '-parin {2}/segreg.txt'
+                  .format(template_dir, tmps_dir, param_dir, fr))
+        
+        os.system('msnreg 2 '
+                  '{0}/LVendo_{3}.vtk '
+                  '{0}/LVepi_{3}.vtk '
+                  '{1}/N_LVendo_{3}.vtk '
+                  '{1}/N_LVepi_{3}.vtk '
+                  '-dofin {2}/lv_{3}_nreg.dof.gz '
+                  '-dofout {2}/lv{3}final.dof.gz '
+                  '-ds 4 -symmetric'
+                  .format(template_dir, vtks_dir, tmps_dir, fr))
+    
+        # same number of points    
+        os.system('cardiacsurfacemap '
+                  '{0}/LVendo_{3}.vtk '
+                  '{1}/N_LVendo_{3}.vtk '
+                  '{2}/lv{3}final.dof.gz '
+                  '{1}/F_LVendo_{3}.vtk'
+                  .format(template_dir, vtks_dir, tmps_dir, fr))
+    
+        os.system('cardiacsurfacemap '
+                  '{0}/LVepi_{3}.vtk '
+                  '{1}/N_LVepi_{3}.vtk '
+                  '{2}/lv{3}final.dof.gz '
+                  '{1}/F_LVepi_{3}.vtk'
+                  .format(template_dir, vtks_dir, tmps_dir, fr))
+    
+        os.system('ptransformation '
+                  '{0}/LVmyo_{3}.vtk '
+                  '{1}/F_LVmyo_{3}.vtk '
+                  '-dofin {2}/lv{3}final.dof.gz'
+                  .format(template_dir, vtks_dir, tmps_dir, fr))
+        
+        os.system('cardiacsurfacemap '
+                  '{0}/RV_{3}.vtk '
+                  '{1}/N_RV_{3}.vtk '
+                  '{2}/rv{3}ds8.dof.gz '
+                  '{1}/C_RV_{3}.vtk'
+                  .format(template_dir, vtks_dir, tmps_dir, fr))
+     
+        ###########################################################################
+        os.system('cp {0}/F_LVendo_{1}.vtk {0}/S_LVendo_{1}.vtk'.format(vtks_dir, fr))
+        os.system('cp {0}/F_LVepi_{1}.vtk {0}/S_LVepi_{1}.vtk'.format(vtks_dir, fr))
+        os.system('cp {0}/F_LVmyo_{1}.vtk {0}/S_LVmyo_{1}.vtk'.format(vtks_dir, fr))
+        os.system('cp {0}/F_LVmyo_{1}.vtk {0}/C_LVmyo_{1}.vtk'.format(vtks_dir, fr))
+        os.system('cp {0}/F_LVmyo_{1}.vtk {0}/W_LVmyo_{1}.vtk'.format(vtks_dir, fr))
+        os.system('cp {0}/C_RV_{1}.vtk {0}/S_RV_{1}.vtk'.format(vtks_dir, fr))
+        os.system('cp {0}/C_RV_{1}.vtk {0}/W_RV_{1}.vtk'.format(vtks_dir, fr))
+        
+        ###########################################################################
+        # compute the quantities of the heart with respect to template
+        os.system('cardiacwallthickness '
+                  '{0}/F_LVendo_{1}.vtk '
+                  '{0}/F_LVepi_{1}.vtk '
+                  '-myocardium '
+                  '{0}/W_LVmyo_{1}.vtk'
+                  .format(vtks_dir, fr))
+    
+        os.system('cardiacenlargedistance '
+                  '{0}/S_LVendo_{2}.vtk '
+                  '{0}/S_LVepi_{2}.vtk '
+                  '{1}/LVendo_{2}.vtk '
+                  '{1}/LVepi_{2}.vtk '
+                  '-myocardium '
+                  '{0}/S_LVmyo_{2}.vtk'
+                  .format(vtks_dir, template_dir, fr))
+        
+        os.system('DiscreteCurvatureEstimator '
+                  '{0}/C_LVmyo_{1}.vtk '
+                  '{0}/FC_LVmyo_{1}.vtk'
+                  .format(vtks_dir, fr))
+        
+        os.system('cardiaccurvature '
+                  '{0}/FC_LVmyo_{1}.vtk '
+                  '{0}/C_LVmyo_{1}.vtk '
+                  '-smooth 64'
+                  .format(vtks_dir, fr))
+        
+        os.system('DiscreteCurvatureEstimator '
+                  '{0}/C_RV_{1}.vtk '
+                  '{0}/FC_RV_{1}.vtk'
+                  .format(vtks_dir, fr))
+        
+        os.system('cardiaccurvature '
+                  '{0}/FC_RV_{1}.vtk '
+                  '{0}/C_RV_{1}.vtk '
+                  '-smooth 64'
+                  .format(vtks_dir, fr))
+        
+        os.system('sevaluation '
+                  '{0}/S_RV_{2}.vtk '
+                  '{1}/RV_{2}.vtk '
+                  '-scalar '
+                  '-signed'
+                  .format(vtks_dir, template_dir, fr))
+        
+        os.system('cardiacwallthickness '
+                  '{0}/W_RV_{1}.vtk '
+                  '{0}/N_RVepi_{1}.vtk'
+                  .format(vtks_dir, fr))
+        
+        ###########################################################################
+#        os.system('vtk2txt {0}/C_RV_{1}.vtk {2}/rv_{1}_curvature.txt'.format(vtks_dir, fr, subject_dir))
+#        os.system('vtk2txt {0}/W_RV_{1}.vtk {2}/rv_{1}_wallthickness.txt'.format(vtks_dir, fr, subject_dir))
+#        os.system('vtk2txt {0}/S_RV_{1}.vtk {2}/rv_{1}_signeddistances.txt'.format(vtks_dir, fr, subject_dir))
+#        os.system('vtk2txt {0}/W_LVmyo_{1}.vtk {2}/lv_myo{1}_wallthickness.txt'.format(vtks_dir, fr, subject_dir))
+#        os.system('vtk2txt {0}/C_LVmyo_{1}.vtk {2}/lv_myo{1}_curvature.txt'.format(vtks_dir, fr, subject_dir))
+#        os.system('vtk2txt {0}/S_LVmyo_{1}.vtk {2}/lv_myo{1}_signeddistances.txt'.format(vtks_dir, fr, subject_dir))
+        if fr == 'ED':
+            fr_ = 'ed'
+            os.system('vtk2txt {0}/C_RV_{1}.vtk {2}/rv_{3}_curvature.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/W_RV_{1}.vtk {2}/rv_{3}_wallthickness.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/S_RV_{1}.vtk {2}/rv_{3}_signeddistances.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/W_LVmyo_{1}.vtk {2}/lv_myo{3}_wallthickness.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/C_LVmyo_{1}.vtk {2}/lv_myo{3}_curvature.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/S_LVmyo_{1}.vtk {2}/lv_myo{3}_signeddistances.txt'.format(vtks_dir, fr, subject_dir, fr_))
+        if fr == 'ES':
+            fr_ = 'es'
+            os.system('vtk2txt {0}/C_RV_{1}.vtk {2}/rv_{3}_curvature.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/W_RV_{1}.vtk {2}/rv_{3}_wallthickness.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/S_RV_{1}.vtk {2}/rv_{3}_signeddistances.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/W_LVmyo_{1}.vtk {2}/lv_myo{3}_wallthickness.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/C_LVmyo_{1}.vtk {2}/lv_myo{3}_curvature.txt'.format(vtks_dir, fr, subject_dir, fr_))
+            os.system('vtk2txt {0}/S_LVmyo_{1}.vtk {2}/lv_myo{3}_signeddistances.txt'.format(vtks_dir, fr, subject_dir, fr_))
+  
+      
+def apply_PC(subject, data_dir, param_dir, template_dir, mirtk):
+    
+    print('  co-registering {0}'.format(subject))    
+    
+    subject_dir = os.path.join(data_dir, subject)
+        
+    if os.path.isdir(subject_dir):
+        
+        tmps_dir = '{0}/tmps'.format(subject_dir)
+
+        vtks_dir = '{0}/vtks'.format(subject_dir)
+
+        dofs_dir = '{0}/dofs'.format(subject_dir)
+
+        meshGeneration(subject_dir, template_dir, param_dir, tmps_dir, vtks_dir, dofs_dir)
+        
+        print('  finish generating meshes from segmentations in {0}'.format(subject))
+        
+    else:  
+        print('  {0} is not a valid directory, do nothing'.format(subject_dir))
+
+
+def meshCoregstration(dir_0, dir_2, dir_3, coreNo, parallel, mirtk):
+               
+    if parallel:
+    
+        print('Generate meshes from segmentations running on {0} cores'.format(coreNo))
+        
+        pool = Pool(processes = coreNo) 
+        
+        # partial only in Python 2.7+
+        pool.map(partial(apply_PC, 
+                         data_dir=dir_0,  
+                         param_dir=dir_2, 
+                         template_dir=dir_3,
+                         mirtk=mirtk), 
+                         sorted(os.listdir(dir_0)))       
+                
+    else:
+        
+        print('Generate meshes from segmentations running subsequently')
+                
+        data_dir, param_dir, template_dir = dir_0, dir_2, dir_3
+                                     
+        for subject in sorted(os.listdir(data_dir)):
+            
+            print('  co-registering {0}'.format(subject))    
+            
+            subject_dir = os.path.join(data_dir, subject)
+
+            if not os.path.isdir(subject_dir):
+                
+                print('  {0} is not a valid folder, Skip'.format(subject_dir))
+                
+                continue 
+                   
+            tmps_dir = '{0}/tmps'.format(subject_dir)
+
+            vtks_dir = '{0}/vtks'.format(subject_dir)
+
+            dofs_dir = '{0}/dofs'.format(subject_dir)
+
+            meshGeneration(subject_dir, template_dir, param_dir, tmps_dir, vtks_dir, dofs_dir)
+            
+            print('  finish generating meshes from segmentations in {0}'.format(subject))
\ No newline at end of file