Merge branch 'master' of https://gitlab.pks.mpg.de/mesoscopic-physics-of-life/frap_theory

Auto_mesh/Fast_auto_mesh.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import tifffile\n",
+    "\n",
+    "from frap_img import find_sigma, contour_from_tif\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "# find_sigma(tif[9][20], tif[10][20],0.23) # we need to play with sigma until we find a value at which all the contours are black and the rest is not\n",
+    "\n",
+    "\n",
+    "\n",
+    "base_folder = 'Analysis/';\n",
+    "tif = tifffile.imread(base_folder + 'halfFrapBig_t000000_cropped.tif')\n",
+    "fig = plt.figure()\n",
+    "x_mid, y_mid, z_mid, r = np.genfromtxt(base_folder + '/mid_point_and_radius.csv', delimiter=',')\n",
+    "\n",
+    "z = 0.3\n",
+    "number_pixel = 207\n",
+    "sigma = 0.23\n",
+    "d = 20\n",
+    "FRAP = 1\n",
+    "N = 8 # when FRAP lot of point in the midlle are created => can lead to a long time for meshing 1D. The number of points is divided by this value\n",
+    "\n",
+    "# Recommended values of precision = [best:0.0025, worst:0.05]\n",
+    "lc_min = 0.032\n",
+    "lc_min_selected = 0.01\n",
+    "lc_min_frap = 0.02\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "contour_from_tif(FRAP,tif[9][20], tif[10][20],sigma,d,'workfilefast.geo', lc_min, lc_min_selected, lc_min_frap,z,number_pixel,N)\n",
+    "\n",
+    "N = 14\n",
+    "# contour_from_tif(FRAP,tif[9][30], tif[10][30],sigma,d,'workfile.geo', lc_min, lc_min_selected, lc_min_frap,z,number_pixel,N)\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# "
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

Auto_mesh/frap_img.py

+from collections import OrderedDict
+from copy import deepcopy
+from math import *
+from scipy import ndimage, misc
+from scipy.interpolate import RegularGridInterpolator
+from scipy.signal import savgol_filter
+from scipy.spatial import distance
+import matplotlib.pyplot as plt
+import numpy as np
+import tifffile
+
+
+
+
+def find_sigma(tif_frame, tif_frame_after_frap ,sigma):
+
+
+    result = ndimage.gaussian_laplace(tif_frame,sigma)
+    plt.imshow(tif_frame)
+    plt.title('Experimental data')
+    bar = plt.colorbar() 
+    bar.set_label('Intensity',size=14)
+    plt.show()
+    plt.imshow(result)
+    plt.title('Gaussian Filter and Laplace filter')
+    bar = plt.colorbar() 
+    bar.set_label('Intensity',size=14)
+    plt.show()
+    
+    
+    result2 = ndimage.gaussian_laplace(tif_frame_after_frap,sigma)
+    plt.imshow(tif_frame_after_frap)
+    plt.title('Experimental data')
+    bar = plt.colorbar() 
+    bar.set_label('Intensity',size=14)
+    plt.show()
+    plt.imshow(result2)
+    plt.title('Gaussian Filter and Laplace filter')
+    bar = plt.colorbar() 
+    bar.set_label('Intensity',size=14)
+    plt.show()
+
+
+
+
+
+def contour_from_tif(FRAP,tif_frame, tif_frame_after_frap ,sigma, d, mesh_file, lc_min, lc_min_selected, lc_min_frap,z,number_pixel,N):
+        
+    # mesh parameters all contours
+    lc_max = 1
+    d_min_mesh = 0.03
+    d_max_mesh = 0.3
+    
+    # mesh parameters selected contours
+    lc_max_selected = 0.5
+    d_min_mesh_selected = 0.002
+    d_max_mesh_selected = 1
+    
+    # mesh parameters FRAP contours
+    lc_max_frap = 0.5
+    d_min_frap = 0.002
+    d_max_frap = 1
+    
+    result = ndimage.gaussian_laplace(tif_frame,sigma)
+    plt.imshow(tif_frame)
+    plt.title('Experimental data')
+    bar = plt.colorbar() 
+    bar.set_label('Intensity',size=14)
+    plt.show()
+    plt.imshow(result)
+    plt.title('Gaussian Filter and Laplace filter')
+    bar = plt.colorbar() 
+    bar.set_label('Intensity',size=14)
+    plt.show()
+    
+    if FRAP == 1:
+    
+        result2 = ndimage.gaussian_laplace(tif_frame_after_frap,sigma)
+        plt.imshow(tif_frame_after_frap)
+        plt.title('Experimental data')
+        bar = plt.colorbar() 
+        bar.set_label('Intensity',size=14)
+        plt.show()
+        plt.imshow(result2)
+        plt.title('Gaussian Filter and Laplace filter')
+        bar = plt.colorbar() 
+        bar.set_label('Intensity',size=14)
+        plt.show()
+    
+    
+    MIN =np.min(result)
+    Mean = np.mean(result)
+
+    Contours_x = []
+    Contours_y = []
+    
+    # thresold
+    result[result > (Mean-MIN)/2] = 0
+    Contours_y, Contours_x = np.nonzero(result)
+
+    
+    coords = np.eye(int(len(Contours_x)),2)
+    
+    for i in range(int(len(Contours_x))):
+        
+        coords[i][0] = Contours_x[i]
+        coords[i][1] = Contours_y[i]
+        
+    coords_bis = deepcopy(coords)
+    Poly_x = []
+    LL = []
+    while coords.shape[0] != 0:
+        count = 0
+        p_x = []
+        sP = coords
+        pA = coords[0,:] 
+        distances = np.linalg.norm(sP - pA, ord=2, axis=1.)  # 'distances' is a list
+        LL = [i for i,x in enumerate(distances) if x<=d]
+        GG = [coords[i] for i in LL]        
+        if len(GG) != 0:
+            A = ((coords_bis[:,None] == GG).all(2).any(1))
+            result = np.where(A == True)
+            p_x = p_x + result[0].tolist()
+            p_x = list(OrderedDict.fromkeys(p_x))
+            coords = np.delete(coords, LL, axis=0)
+        while count != len(p_x) :
+            c = 0
+            sP = coords
+            if coords.shape[0] != 0:
+                pA = coords_bis[p_x[count],:]                
+                distances = np.linalg.norm(sP - pA, ord=2, axis=1.)  # 'distances' is a list
+                LL = [i for i,x in enumerate(distances) if x<=d]
+                GG = [coords[i] for i in LL]
+                if len(GG) != 0:
+                    A = ((coords_bis[:,None] == GG).all(2).any(1))
+                    A = np.array(A)
+                    result = np.where(A == True)
+                    p_x = p_x + result[0].tolist()
+                    p_x = list(OrderedDict.fromkeys(p_x))
+                    coords = np.delete(coords, LL, axis=0)
+            count = count +1
+        Poly_x.append(p_x)
+
+    Y = []
+    X = []
+    for j in range(len(Poly_x)):
+        for i in range(len(Poly_x[j])):
+            Y.append(coords_bis[Poly_x[j][i]][0])
+            X.append(coords_bis[Poly_x[j][i]][1])
+        plt.scatter(Y,X, label = "%i" %j)
+
+        Y = []
+        X = []
+
+    plt.axis([0,207, 207, 0])
+    plt.gca().set_aspect('equal', adjustable='box')
+    plt.legend()
+    plt.title('Contours')
+    plt.show()
+        
+    print('Choose a contours:') 
+    choose = input()     
+    
+    if FRAP == 1:
+        MIN =np.min(result2)
+        Mean = np.mean(result2)
+
+        Contours_x = []
+        Contours_y = []
+    
+        # thresold
+        result2[result2 > (Mean-MIN)/2] = 0
+        Contours_y, Contours_x = np.nonzero(result2)
+
+    
+        coords2 = np.eye(int(len(Contours_x)),2)    
+        for i in range(int(len(Contours_x))):
+        
+            coords2[i][0] = Contours_x[i]
+            coords2[i][1] = Contours_y[i]
+        
+        coords_bis2 = deepcopy(coords2)
+        
+        
+        Poly_x2 = []           
+        LL = []
+        count = 0
+
+        
+        p_x = []
+        sP = coords2
+        pA = coords_bis[Poly_x[int(choose)][0]][:]
+        distances = np.linalg.norm(sP - pA, ord=2, axis=1.)  # 'distances' is a list
+        LL = [i for i,x in enumerate(distances) if x<=d]
+        GG = [coords2[i] for i in LL]        
+        if len(GG) != 0:
+            A = ((coords_bis2[:,None] == GG).all(2).any(1))
+            A = np.array(A)
+            result = np.where(A == True)
+            p_x = p_x + result[0].tolist()
+            p_x = list(OrderedDict.fromkeys(p_x))
+            coords2 = np.delete(coords2, LL, axis=0)
+        while count != len(p_x) :
+            c = 0
+            sP = coords2
+            if coords2.shape[0] != 0:
+                pA = coords_bis2[p_x[count],:]                
+                distances = np.linalg.norm(sP - pA, ord=2, axis=1.)  # 'distances' is a list
+                LL = [i for i,x in enumerate(distances) if x<=d]
+                GG = [coords2[i] for i in LL]
+                if len(GG) != 0:
+                    A = ((coords_bis2[:,None] == GG).all(2).any(1))
+                    A = np.array(A)
+                    result = np.where(A == True)
+                    p_x = p_x + result[0].tolist()
+                    p_x = list(OrderedDict.fromkeys(p_x))
+                    coords2 = np.delete(coords2, LL, axis=0)
+            count = count +1
+        Y = []
+        X = []
+        for j in range(len(Poly_x)):
+            if j == int(choose):
+                for k in range(0, int(len(p_x)-1-N),N):
+                    Y.append(coords_bis2[p_x[k]][0])
+                    X.append(coords_bis2[p_x[k]][1])
+            for i in range(len(Poly_x[j])):
+                Y.append(coords_bis[Poly_x[j][i]][0])
+                X.append(coords_bis[Poly_x[j][i]][1])
+            plt.scatter(Y,X, label = "%i" %j)
+            Y = []
+            X = []
+        plt.axis([0,207, 207, 0])
+        plt.gca().set_aspect('equal', adjustable='box')
+        plt.legend()
+        plt.title('Contours')
+        plt.show()
+
+    # We create the file
+    f = open(mesh_file, 'w')
+    # Creation of the geometry
+    f.write('//+\n' )
+    f.write('lc = 1;\n' )
+    f.write('//+\n' )
+    f.write('Point(4) = {0, 1, 0, lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(3) = {1, 1, 0, lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(2) = {1, 0, 0, lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(1) = {0, 0, 0, lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(8) = {0, 1, '+ str(z)+', lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(7) = {1, 1, '+ str(z)+', lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(6) = {1, 0, '+ str(z)+', lc};\n' )
+    f.write('//+\n' )
+    f.write('Point(5) = {0, 0, '+ str(z)+', lc};\n' )
+    f.write('//+\n' )
+    f.write('Line(6) = {2, 1};\n' )
+    f.write('//+\n' )
+    f.write('Line(5) = {3, 2};\n' )
+    f.write('//+\n' )
+    f.write('Line(8) = {4, 3};\n' )
+    f.write('//+\n' )
+    f.write('Line(7) = {1, 4};\n' )
+    f.write('//+\n' )
+    f.write('Line(3) = {6, 5};\n' )
+    f.write('//+\n' )
+    f.write('Line(2) = {7, 6};\n' )
+    f.write('//+\n' )
+    f.write('Line(1) = {8, 7};\n' )
+    f.write('//+\n' )
+    f.write('Line(4) = {5, 8};\n' )
+    f.write('//+\n' )
+    f.write('Line(30000000) = {5, 1};\n' )
+    f.write('//+\n' )
+    f.write('Line(10000000) = {2, 6};\n' )
+    f.write('//+\n' )
+    f.write('Line(9) = {3, 7};\n' )
+    f.write('//+\n' )
+    f.write('Line(20000000) = {8, 4};\n' )
+    f.write('//+\n' )
+    f.write('Line Loop(13) = {9, 2, -10000000, -5};\n' )
+    f.write('//+\n' )
+    f.write('Plane Surface(14) = {13};\n' )
+    f.write('//+\n' )
+    f.write('Line Loop(15) = {1, -9, -8, -20000000};\n' )
+    f.write('//+\n' )
+    f.write('Plane Surface(16) = {15};\n' )
+    f.write('//+\n' )
+    f.write('Line Loop(17) = {8, 5, 6, 7};\n' )
+    f.write('//+\n' )
+    f.write('Plane Surface(18) = {17};\n' )
+    f.write('//+\n' )
+    f.write('Line Loop(19) = {3, 30000000, -6, 10000000};\n' )
+    f.write('//+\n' )
+    f.write('Plane Surface(20) = {19};\n' )
+    f.write('//+\n' )
+    f.write('Line Loop(21) = {30000000, 7, -20000000, -4};\n' )
+    f.write('//+\n' )
+    f.write('Plane Surface(22) = {21};\n' )
+    f.write('//+\n' )
+    f.write('Line Loop(23) ={2, 3, 4, 1};\n' )
+    f.write('//+\n' )
+    f.write('Plane Surface(24) = {-23};\n' )
+    f.write('//+\n' )
+    f.write('Surface Loop(25) = {24, 14, 16, 18, 20, 22};\n' )
+    f.write('//+\n' )
+    f.write('Volume(26) = {25};\n' )
+    f.write('//+\n' )
+    f.write('Physical Surface(27) = {16};\n' )
+    f.write('//+\n' )
+    f.write('Physical Surface(28) = {20};\n' )
+    f.write('//+\n' )
+    f.write('Physical Volume(29) = {26};\n' )
+    f.write('//+\n' )
+    f.write('Physical Surface(30) = {24};\n' )
+    f.write('//+\n' )
+    f.write('Physical Surface(31) = {18};\n' )
+    f.write('//+\n' )
+    
+    # Creation of the FRAP points
+    if FRAP ==1:
+        for j in range(0,int((len(p_x))),N):
+
+            f.write('Point(3000'+str(j)+') = {'+str(((coords_bis2[p_x[j]][0]/number_pixel)))+', '+ str((1-(coords_bis2[p_x[j]][1]/number_pixel)))+', 0, lc};\n' )
+            f.write('//+\n' )
+
+            f.write('Point('+str(4000*int(z+1))+str(j)+') = {'+str(((coords_bis2[p_x[j]][0]/number_pixel)))+', '+ str((1-(coords_bis2[p_x[j]][1]/number_pixel)))+','+ str(z)+', lc};\n' )
+            f.write('//+\n' )
+            
+    # Creation of the contours points
+    for i in range(len(Poly_x)):
+        for j in range(len(Poly_x[i])):
+
+            f.write('Point('+str(i+1)+str(j)+') = {'+str(((coords_bis[Poly_x[i][j]][0]/number_pixel)))+', '+ str((1-(coords_bis[Poly_x[i][j]][1]/number_pixel)))+', 0, lc};\n' )
+            f.write('//+\n' )
+            
+            f.write('Point('+str(1000*int(z+1))+str(i+1)+str(j)+') = {'+str(((coords_bis[Poly_x[i][j]][0]/number_pixel)))+', '+ str((1-(coords_bis[Poly_x[i][j]][1]/number_pixel)))+','+ str(z)+', lc};\n' )
+            f.write('//+\n' )
+            
+    # Creation of the contours Lines
+    for i in range(len(Poly_x)):
+        for j in range(len(Poly_x[i])-1):
+            if sqrt((coords_bis[Poly_x[i][j]][1] - coords_bis[Poly_x[i][j+1]][1])*(coords_bis[Poly_x[i][j]][1] -coords_bis[Poly_x[i][j+1]][1]) + ((coords_bis[Poly_x[i][j]][0] -  coords_bis[Poly_x[i][j+1]][0])*(coords_bis[Poly_x[i][j]][0] - coords_bis[Poly_x[i][j+1]][0])))< d :
+                f.write('Line('+str(i+1)+str(j)+') = {'+str(i+1)+str(j)+','+str(i+1)+str(j+1)+'};\n' )
+                f.write('//+\n' )
+                
+                f.write('Line('+str(1000*int(z+1))+str(i+1)+str(j)+') = {'+str(1000*int(z+1))+str(i+1)+str(j)+','+str(1000*int(z+1))+str(i+1)+str(j+1)+'};\n' )
+                f.write('//+\n' )
+                
+                f.write('Line('+str(2000*int(z+1))+str(i+1)+str(j)+') = {'+str(1000*int(z+1))+str(i+1)+str(j)+','+str(i+1)+str(j)+'};\n' )
+                f.write('//+\n' )
+                
+    # Creation of the FRAP Lines
+    
+    
+    if FRAP ==1:
+        for i in range(0,int(len(p_x)-N-1),N):
+            if sqrt((coords_bis2[p_x[i]][1] - coords_bis2[p_x[i+N]][1])*(coords_bis2[p_x[i]][1] - coords_bis2[p_x[i+N]][1]) + ((coords_bis2[p_x[i]][0] -  coords_bis2[p_x[i+N]][0])*(coords_bis2[p_x[i]][0] - coords_bis2[p_x[i+N]][0])))< d :
+
+                f.write('Line(3000'+str(i+1)+') = {3000'+str(i)+','+'3000'+str(i+N)+'};\n' )
+                f.write('//+\n' )
+
+                f.write('Line(4000'+str(i+1)+') = {'+str(4000*int(z+1))+str(i)+','+str(4000*int(z+1))+str(i+N)+'};\n' )
+                f.write('//+\n' )
+
+
+                f.write('Line(5000'+str(i+1)+') = {'+str(4000*int(z+1))+str(i)+','+'3000'+str(i+N)+'};\n' )
+                f.write('//+\n' )
+                
+    # Creation of the contours mesh
+    f.write('Field[1] = Distance;\n' )
+    f.write('//+\n' )
+    f.write('Field[1].EdgesList = {' )
+    
+    
+    
+          
+    
+    
+    
+    
+    for i in range(len(Poly_x)):
+        for j in range(len(Poly_x[i])-1):
+            if sqrt((coords_bis[Poly_x[i][j]][1] - coords_bis[Poly_x[i][j+1]][1])*(coords_bis[Poly_x[i][j]][1] -coords_bis[Poly_x[i][j+1]][1]) + ((coords_bis[Poly_x[i][j]][0] -  coords_bis[Poly_x[i][j+1]][0])*(coords_bis[Poly_x[i][j]][0] - coords_bis[Poly_x[i][j+1]][0])))< d :
+                
+                f.write(str(1000*int(z+1))+str(i+1)+str(j)+',')
+                f.write(str(2000*int(z+1))+str(i+1)+str(j)+',')                
+                f.write(str(i+1)+str(j)+',')
+                a = i+1
+                b= j
+    f.write(str(a)+str(b))
+    f.write('};\n' )
+    f.write('//+\n' )
+    f.write('Field[1].NNodesByEdge = 100;\n' )
+    f.write('//+\n' )
+    f.write('Field[2] = Threshold;\n' )
+    f.write('//+\n' )
+    f.write('Field[2].IField = 1;\n' )
+    f.write('//+\n' )
+    f.write('Field[2].LcMin = '+ str(lc_min)+';\n' )
+    f.write('//+\n' )
+    f.write('Field[2].LcMax = '+ str(lc_max)+';\n' )
+    f.write('//+\n' )
+    f.write('Field[2].DistMin = '+ str(d_min_mesh)+';\n' )
+    f.write('//+\n' )
+    f.write('Field[2].DistMax = '+ str(d_max_mesh)+';\n' )
+    f.write('//+\n' ) 
+    
+    # Creation of the selected contours mesh
+    f.write('Field[3] = Distance;\n' )
+    f.write('//+\n' )
+    f.write('Field[3].EdgesList = {' )
+    
+    for j in range(len(Poly_x[int(choose)])-1):
+        if sqrt((coords_bis[Poly_x[int(choose)][j]][1] - coords_bis[Poly_x[int(choose)][j+1]][1])*(coords_bis[Poly_x[int(choose)][j]][1] -coords_bis[Poly_x[int(choose)][j+1]][1]) + (coords_bis[Poly_x[int(choose)][j]][0] -  coords_bis[Poly_x[int(choose)][j+1]][0])*(coords_bis[Poly_x[int(choose)][j]][0] - coords_bis[Poly_x[int(choose)][j+1]][0]))< d :
+            f.write(str(int(choose)+1)+str(j)+',')
+            f.write(str(2000*int(z+1))+str(int(choose)+1)+str(j)+',')                
+            f.write(str(1000*int(z+1))+str(int(choose)+1)+str(j)+',')
+            a = int(choose)+1
+            b= j
+    f.write(str(a)+str(b))
+    f.write('};\n' )
+    f.write('//+\n' )
+    f.write('Field[3].NNodesByEdge = 100;\n' )
+    f.write('//+\n' )
+    f.write('Field[4] = Threshold;\n' )
+    f.write('//+\n' )
+    f.write('Field[4].IField = 3;\n' )
+    f.write('//+\n' )
+    f.write('Field[4].LcMin = '+ str(lc_min_selected)+';\n' )
+    f.write('//+\n' )
+    f.write('Field[4].LcMax = '+ str(lc_max_selected)+';\n' )
+    f.write('//+\n' )
+    f.write('Field[4].DistMin = '+ str(d_min_mesh_selected)+';\n' )
+    f.write('//+\n' )
+    f.write('Field[4].DistMax = '+ str(d_max_mesh_selected)+';\n' )
+    f.write('//+\n' )
+    
+    # Creation of the FRAP contours mesh
+    if FRAP == 1:
+        f.write('Field[5] = Distance;\n' )
+        f.write('//+\n' )
+        f.write('Field[5].EdgesList = {' )
+        for i in range(0,int(len(p_x)-N-1),N):
+            f.write('3000'+str(i+1)+',')
+            f.write('4000'+str(i+1)+',')                
+            f.write('5000'+str(i+1)+',')
+            a = 1
+            b= i
+        f.write(str(a)+str(b))
+        f.write('};\n' )
+        f.write('//+\n' )
+        f.write('Field[5].NNodesByEdge = 100;\n' )
+        f.write('//+\n' )
+        f.write('Field[6] = Threshold;\n' )
+        f.write('//+\n' )
+        f.write('Field[6].IField = 5;\n' )
+        f.write('//+\n' )
+        f.write('Field[6].LcMin = '+ str(lc_min_frap)+';\n' )
+        f.write('//+\n' )
+        f.write('Field[6].LcMax = '+ str(lc_max_frap)+';\n' )
+        f.write('//+\n' )
+        f.write('Field[6].DistMin = '+ str(d_min_frap)+';\n' )
+        f.write('//+\n' )
+        f.write('Field[6].DistMax = '+ str(d_max_frap)+';\n' )
+        f.write('//+\n' )
+    f.write('Field[7] = Min;\n' )
+    f.write('//+\n' )
+    if FRAP ==0:
+        f.write('Field[7].FieldsList = {4,2};\n' )
+        f.write('//+\n' )
+        f.write('Background Field = 7;\n' )
+        f.write('//+\n' )
+        f.close()
+    if FRAP ==1:
+        f.write('Field[7].FieldsList = {4,2,6};\n' )
+        f.write('//+\n' )
+        f.write('Background Field = 7;\n' )
+        f.write('//+\n' )
+        f.close()
+#     f = open('workfile.geo', 'r')
+#     f.read()
+#     for line in f:
+#         print(line, end='')
+#     f.close()
+    
+