--- a
+++ b/utils/myocardial_strain.py
@@ -0,0 +1,654 @@
+import numpy as np
+from scipy.ndimage import rotate
+from scipy.ndimage import gaussian_filter
+from scipy.interpolate import interp1d, interp2d
+from scipy.ndimage.measurements import center_of_mass
+
+class Tensor():
+    
+    def __init__(self, Exx, Exy, Exz, Eyx, Eyy, Eyz, Ezx, Ezy, Ezz):
+        
+        self.E1, self.E2, self.E3 = Exx.copy(), Exy.copy(), Exz.copy()
+        self.E4, self.E5, self.E6 = Eyx.copy(), Eyy.copy(), Eyz.copy()
+        self.E7, self.E8, self.E9 = Ezx.copy(), Ezy.copy(), Ezz.copy()
+
+    def asmat(self):
+        return np.array([[self.E1,self.E2,self.E3],
+                         [self.E4,self.E5,self.E6],
+                         [self.E7,self.E8,self.E9]]).transpose((2,3,4,0,1))
+        
+    def asvoigt(self):
+        return self.E1, self.E2, self.E3, self.E4, self.E5, self.E6, self.E7, self.E8, self.E9 
+        
+    def transpose(self):
+        return Tensor(self.E1, self.E4, self.E7, self.E2, self.E5, self.E8, self.E3, self.E6, self.E9)
+    
+    def identity_add(self):
+        self.E1 += 1; self.E5 += 1; self.E9 += 1 
+    
+    def identity_subtract(self):
+        self.E1 -= 1; self.E5 -= 1; self.E9 -= 1 
+        
+    @staticmethod
+    def dot(X, Y):
+        
+        X1,X2,X3,X4,X5,X6,X7,X8,X9=X.asvoigt()
+        Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8,Y9=Y.asvoigt()
+        
+        Z1, Z2, Z3 = X1*Y1 + X2*Y4 + X3*Y7, X1*Y2 + X2*Y5 + X3*Y8, X1*Y3 + X2*Y6 + X3*Y9
+        Z4, Z5, Z6 = X4*Y1 + X5*Y4 + X6*Y7, X4*Y2 + X5*Y5 + X6*Y8, X4*Y3 + X5*Y6 + X6*Y9
+        Z7, Z8, Z9 = X7*Y1 + X8*Y4 + X9*Y7, X7*Y2 + X8*Y5 + X9*Y8, X7*Y3 + X8*Y6 + X9*Y9
+        
+        return Tensor(Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8,Z9)
+    
+
+class MyocardialStrain():
+    
+    def __init__(self, mask, flow):
+                
+        self.mask  = mask
+        self.flow  = flow
+        
+        assert len(mask.shape) == 3
+        assert len(flow.shape) == 4
+        assert mask.shape == flow.shape[:3]
+        assert flow.shape[-1] == 3
+
+    def calculate_strain(self, lv_label=3):
+        
+        cx, cy, cz = center_of_mass(self.mask==lv_label)
+        nx, ny, nz = self.mask.shape
+        
+        self.flow_rot = roll_to_center(self.flow, cx, cy)
+        self.mask_rot = roll_to_center(self.mask, cx, cy)
+
+        ux, uy, uz  = np.array_split(self.flow_rot, 3, -1)
+        Uxx,Uxy,Uxz = np.gradient(np.squeeze(ux))
+        Uyx,Uyy,Uyz = np.gradient(np.squeeze(uy))
+        Uzx,Uzy,Uzz = np.gradient(np.squeeze(uz))
+        
+        F = Tensor(Uxx,Uxy,Uxz,Uyx,Uyy,Uyz,Uzx,Uzy,Uzz)
+
+        F.identity_add()
+        F = F.dot(F.transpose(), F) 
+        F.identity_subtract()
+        
+        self.Err, self.Ecc, self.Erc, self.Ecr = convert_to_polar(mask=self.mask_rot, E=0.5*F.asmat()[:,:,:,:2,:2])
+        
+
+def roll(x, rx, ry):
+    x = np.roll(x, rx, axis=0)
+    return np.roll(x, ry, axis=1)
+
+def roll_to_center(x, cx, cy):
+    nx, ny = x.shape[:2]
+    return roll(x,  int(nx//2-cx), int(ny//2-cy))
+
+def polar_grid(nx=128, ny=128):
+    x, y = np.meshgrid(np.linspace(-nx//2, nx//2, nx), np.linspace(-ny//2, ny//2, ny))
+    phi  = (np.rad2deg(np.arctan2(y, x)) + 180).T
+    r    = np.sqrt(x**2+y**2+1e-8)
+    return phi, r
+
+def convert_to_polar(mask, E):
+
+    phi = polar_grid(*E.shape[:2])[0]
+    Err = np.zeros(mask.shape)
+    Ecc = np.zeros(mask.shape)
+    Erc = np.zeros(mask.shape)
+    Ecr = np.zeros(mask.shape)
+    for k in range(mask.shape[-1]):
+        cos = np.cos(np.deg2rad(phi))
+        sin = np.sin(np.deg2rad(phi))
+
+        Exx, Exy, Eyx, Eyy = E[:,:,k,0,0],E[:,:,k,0,1],E[:,:,k,1,0],E[:,:,k,1,1]
+        Err[:,:,k] +=  cos*( cos*Exx+sin*Exy) + sin*( cos*Eyx+sin*Eyy)
+        Ecc[:,:,k] += -sin*(-sin*Exx+cos*Exy) + cos*(-sin*Eyx+cos*Eyy)
+        Erc[:,:,k] +=  cos*(-sin*Exx+cos*Exy) + sin*(-sin*Eyx+cos*Eyy)
+        Ecr[:,:,k] += -sin*( cos*Exx+sin*Exy) + cos*( cos*Eyx+sin*Eyy)
+
+    return Err, Ecc, Erc, Ecr
+
+    
+def convert_to_aha4d(tensor, mask):
+    Tensor = tensor.copy()
+    Mask   = mask.copy()
+
+    Tensor[Mask!=2,:] = 0
+
+    # rotate to have RV center of mass on the right
+    angle  = _get_lv2rv_angle(Mask)
+    Tensor = rotate(Tensor,  -angle, reshape=False, order=0)
+    Mask   = rotate(Mask, -angle, reshape=False, order=1).clip(0,3).round()
+
+    # roll to center 
+    cx, cy = center_of_mass(Mask>1)[:2]
+    Tensor = np.flipud(np.rot90(_roll_to_center(Tensor, cx, cy)))
+    Mask   = np.flipud(np.rot90(_roll_to_center(Mask, cx, cy)))
+
+    # remove slices that do not contain tissue labels
+    ID     = Mask.sum(axis=(0,1))>0
+    Mask   = Mask[:,:,ID]
+    Tensor = Tensor[:,:,ID]
+    
+    return Tensor, Mask
+    
+class PolarMap():
+    
+    def __init__(self, Err, Ecc, mask):
+        """ Plot Err and Ecc in a PolarMap using a segmentation of the heart as reference. 
+            Assuming mask==2 yields the myocardium labels. 
+        
+        """
+        
+        self.Err  = Err
+        self.Ecc  = Ecc
+        self.mask = mask
+        
+    def project_to_aha_polar_map(self):
+       
+        Err  = self.Err.copy()
+        Ecc  = self.Ecc.copy()
+        mask = self.mask.copy()
+
+        Err[mask!=2] = 0
+        Ecc[mask!=2] = 0
+
+        # rotate to have RV center of mass on the right
+        angle = _get_lv2rv_angle(mask)
+        Ecc   = rotate(Ecc,  -angle, reshape=False, order=0)
+        Err   = rotate(Err,  -angle, reshape=False, order=0)
+        mask  = rotate(mask, -angle, reshape=False, order=1).clip(0,3).round()
+        
+        # roll to center 
+        cx, cy = center_of_mass(mask>1)[:2]
+        Ecc  = np.flipud(np.rot90(_roll_to_center(Ecc, cx, cy)))
+        Err  = np.flipud(np.rot90(_roll_to_center(Err, cx, cy)))
+        mask = np.flipud(np.rot90(_roll_to_center(mask, cx, cy)))
+        
+        # remove slices that do not contain tissue labels
+        ID   = mask.sum(axis=(0,1))>0
+        mask = mask[:,:,ID]
+        Err  = Err[:,:,ID]
+        Ecc  = Ecc[:,:,ID]
+
+
+        
+
+        Err = Err.transpose((2,0,1))
+        Ecc = Ecc.transpose((2,0,1))
+        print('... radial strain')
+        V_err = self._project_to_aha_polar_map(Err)
+        print('... circumferential strain')
+        V_ecc = self._project_to_aha_polar_map(Ecc)
+
+        results = {'V_err':V_err, 'V_ecc':V_ecc, 'mask':mask}
+
+        return results
+        
+    def _project_to_aha_polar_map(self, E, nphi=360, nrad=100, dphi=1):
+        
+        nz = E.shape[0]
+        angles = np.arange(0, nphi, dphi)
+        V      = np.zeros((nz, angles.shape[0], nrad))
+
+        for rj in range(nz):
+
+            Err_q  = _inpter2(E[rj])
+
+            PHI, R = _polar_grid(*Err_q.shape)
+            PHI = PHI.ravel()
+            R   = R.ravel()
+
+            for k, pmin in enumerate(angles):
+                pmax = pmin + dphi/2.0
+
+                # Get values for angle segment
+                PHI_SEGMENT = (PHI>=pmin)&(PHI<=pmax)
+                Rk   = R[PHI_SEGMENT]
+                PHIk = PHI[PHI_SEGMENT]
+                Vk   = Err_q.ravel()[PHI_SEGMENT]
+                
+                
+                Rk = Rk[np.abs(Vk)!=0]
+                Vk = Vk[np.abs(Vk)!=0]
+                
+                if len(Vk) == 0:
+                    continue # this might not be the best
+                Rk = _rescale_linear(Rk, rj, rj + 1)
+
+                r = np.arange(rj, rj+1, 1.0/nrad)
+                f = interp1d(Rk, Vk)
+                v = f(r)
+
+                V[rj,k] += v
+
+        return V        
+        
+        
+    def construct_polar_map(self, tensor, start=30, stop=70, sigma=12):
+
+        E  = tensor.copy()
+        mu = E[:,:,start:stop].mean()
+
+        nz = E.shape[0]
+        E  = np.concatenate(np.array_split(E[:,:,start:stop], nz), axis=-1)[0]
+
+        old = E.shape[1]/nz*1.
+        for j in range(nz-1):
+            xi = int(old//2+j*old)
+            xj = int(old+old//2+j*old)
+            E[:,xi:xj] = gaussian_filter(E[:,xi:xj],sigma=sigma, mode='wrap')
+            E[:,xi:xj] = gaussian_filter(E[:,xi:xj],sigma=sigma, mode='wrap')
+
+        E = np.stack(np.array_split(E,nz,axis=1))
+
+        E = gaussian_filter(E, sigma=sigma, mode='wrap')
+        E = gaussian_filter(E, sigma=sigma, mode='wrap')
+
+        E = [E[0][None]] + [E[1:3]] + np.array_split(E[3:], 2, axis=0)
+
+        E = [np.mean(E[i], axis=0) for i in range(4)]
+        E = np.concatenate(E, axis=1)
+
+        old = E.shape[1]/4
+        for j in range(3):
+            xi = int(old//2+j*old)
+            xj = int(old+old//2+j*old)
+            E[:,xi:xj] = gaussian_filter(E[:,xi:xj], sigma=sigma, mode='wrap')
+            E[:,xi:xj] = gaussian_filter(E[:,xi:xj], sigma=sigma, mode='wrap')
+
+        E = gaussian_filter(E, sigma=sigma, mode='wrap')
+        E = gaussian_filter(E, sigma=sigma, mode='wrap')
+
+        mu = [mu] + self._get_17segments(E)
+
+        return E, mu 
+    
+    def _get_17segments(self, data):
+        c1,c2,c3,c4 = np.array_split(data,4,axis=-1)
+        c2 = np.roll(c2,-45,0)
+        #c2 = np.roll(c2,-90,0)
+
+        c4 = [np.mean(ci) for ci in np.array_split(c4,6,axis=0)]
+        c4 = list(np.roll(np.array(c4),-1))
+        c3 = [np.mean(ci) for ci in np.array_split(c3,6,axis=0)]
+        c3 = list(np.roll(np.array(c3),-1))
+        c2 = [np.mean(ci) for ci in np.array_split(c2,4,axis=0)]
+        #c2 = list(np.roll(np.array(c2),-1))
+        c1 = [np.mean(c1)]
+
+        c = c4 + c3 + c2 + c1
+
+        return c
+
+    def _get_17segments_RC(self, data1,data2):
+        
+        def _rc(a,b):
+            #return np.mean(np.abs((b-a)/b)*100)
+            return np.mean(((b-a)/b)*100)
+        
+        c1_1,c2_1,c3_1,c4_1 = np.array_split(data1,4,axis=-1)
+        c1_2,c2_2,c3_2,c4_2 = np.array_split(data2,4,axis=-1)
+
+        c4 = [_rc(ci1,ci2) for ci1,ci2 in zip(np.array_split(c4_1,6,axis=0),
+                                              np.array_split(c4_2,6,axis=0))]
+        c4 = list(np.roll(np.array(c4),-1))
+        
+        c3 = [_rc(ci1,ci2) for ci1,ci2 in zip(np.array_split(c3_1,6,axis=0),
+                                              np.array_split(c3_2,6,axis=0))]
+        c3 = list(np.roll(np.array(c3),-1))
+        
+        c2 = [_rc(ci1,ci2) for ci1,ci2 in zip(np.array_split(c2_1,4,axis=0),
+                                              np.array_split(c2_2,4,axis=0))]
+        c2 = list(np.roll(np.array(c2),-1))
+        c1 = [_rc(c1_1,c1_2)]
+
+        c = c4 + c3 + c2 + c1
+
+        return c
+    
+def _roll(x, rx, ry):
+    x = np.roll(x, rx, axis=0)
+    return np.roll(x, ry, axis=1)
+
+def _roll_to_center(x, cx, cy):
+    nx, ny = x.shape[:2]
+    return _roll(x,  int(nx//2-cx), int(ny//2-cy))
+
+def _py_ang(v1, v2):
+    """ Returns the angle in degrees between vectors 'v1' and 'v2'. """
+    cosang = np.dot(v1, v2)
+    sinang = np.linalg.norm(np.cross(v1, v2))
+    return np.rad2deg(np.arctan2(sinang, cosang))
+
+def _polar_grid(nx=128, ny=128):
+    x, y = np.meshgrid(np.linspace(-nx//2, nx//2, nx), np.linspace(-ny//2, ny//2, ny))
+    phi  = (np.rad2deg(np.arctan2(y, x)) + 180).T
+    r    = np.sqrt(x**2+y**2+1e-8)
+    return phi, r
+
+def _rescale_linear(array, new_min, new_max):
+    minimum, maximum = np.min(array), np.max(array)
+    m = (new_max-new_min) / (maximum-minimum)
+    b = new_min - m* minimum
+    return m*array + b
+
+def _inpter2(Eij, k=10):
+    nx, ny = Eij.shape
+    
+    x  = np.linspace(0,nx-1,nx)
+    y  = np.linspace(0,ny-1,ny)
+    xq = np.linspace(0,nx-1,nx*k)
+    yq = np.linspace(0,ny-1,ny*k)
+    
+    f = interp2d(x,y,Eij,kind='linear')
+    
+    return f(xq,yq)    
+    
+def _get_lv2rv_angle(mask):
+    cx_lv, cy_lv = center_of_mass(mask>1)[:2]
+    cx_rv, cy_rv = center_of_mass(mask==1)[:2]
+    phi_angle    = _py_ang([cx_rv-cx_lv, cy_rv-cy_lv], [0, 1])
+    return phi_angle    
+    
+    
+    
+### FUNCTIONS TO PLOT THE POLAR MAP 
+
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+from copy import deepcopy
+
+def _write(ax, mu, j, theta_i, i, width=2):
+    xi, yi = polar2cart(0,  theta_i)
+    xf, yf = polar2cart(35, theta_i)
+
+    l = Line2D([40-xi,40-xf], [40-yi,40-yf], color='black', linewidth=width)
+    ax.add_line(l)
+    xi, yi = polar2cart(30, theta_i + 2*np.pi/12)
+    ax.text(40-xi-.3, 40-yi, '%d' %(mu[j][i]), weight='bold', fontsize=14)
+    
+def write(ax, mu, j, width=2):
+    if j > 1:
+        for i in range(6):
+            theta_i = 2*np.pi - i*60*np.pi/180 + 2*60*np.pi/180
+            _write(ax, mu, j, theta_i, i)
+    if j == 1:
+        for i in range(4):
+            theta_i = i*90*np.pi/180 - 45*np.pi/180
+            _write(ax, mu, j, theta_i, i)
+    if j ==0:
+        ax.text(40-.3, 40, '%d' %(mu[j][0]), weight='bold', fontsize=14)
+
+def plot_bullseye(data,mu,vmin=None,vmax=None, savepath=None,cmap='RdBu_r', label='GPRS (%)', 
+                  std=None,cbar=False,color='white', fs=20, xshift=0, yshift=0, ptype='mesh',frac=False):
+    
+    rho     = np.arange(0,4,4.0/data.shape[1]);
+    Theta   = np.deg2rad(range(data.shape[0]))
+    [th, r] = np.meshgrid(Theta, rho);
+
+    fig, ax = plt.subplots(figsize=(6,6))
+    #fig.subplots_adjust(left=0,right=1,bottom=0,top=1)
+    #ax.axis('tight') creates errors 
+    #ax.axis('off')
+    if ptype == 'mesh':
+        im = ax.pcolormesh(r*np.cos(Theta), r*np.sin(Theta), 100*data.T, 
+                           vmin=vmin,vmax=vmax,cmap=cmap,shading='gouraud')
+    else:
+        im = ax.contourf(r*np.cos(Theta), r*np.sin(Theta), 100*data.T, 
+                           vmin=vmin,vmax=vmax,cmap=cmap,shading='gouraud')
+    if cbar:
+        cbar = plt.colorbar(im, cax=fig.add_axes([0.15, -0.03, 0.7, 0.05]), orientation='horizontal')
+
+        new_ticks = []
+        new_ticks_labels = []
+        for i,tick in enumerate(cbar.ax.get_xticks()):
+            if i % 2 == 0:
+                new_ticks.append(np.round(tick))
+                new_ticks_labels.append(str(int(np.round(tick))))
+
+        cbar.set_ticks(new_ticks);
+        cbar.set_ticklabels(new_ticks_labels);
+
+        # override if vmin is provided, assume vmax is provided too for now
+        if vmin is not None:
+            cbar.set_ticks([vmin, (vmax+vmin)/2.0, vmax]);
+            cbar.set_ticklabels(['%d'%(i) for i in [vmin, (vmax+vmin)/2.0, vmax]]);
+        cbar.ax.tick_params(labelsize=18)
+        cbar.set_label(label, fontsize=26, weight='bold')
+
+    ax.axis('off')
+    if std is not None:
+        draw_circle_group(ax,100*np.array(mu),100*np.array(std))
+    if frac:
+        draw_circle_frac(ax,np.array(mu), color=color)
+    else:
+        draw_circle(ax,100*np.array(mu), color=color, fs=fs, xshift=xshift, yshift=yshift)
+    if savepath is not None:
+        if not cbar:
+            plt.tight_layout()
+        plt.savefig(savepath, dpi=600)
+    plt.show()
+    
+def plot_bullseye_ratio(data,mu,vmin=None,vmax=None, savepath=None,cmap='RdBu_r', label='GPRS (%)', 
+                  std=None,cbar=False,color='white',ptype='mesh',frac=False):
+    
+    rho     = np.arange(0,4,4.0/data.shape[1]);
+    Theta   = np.deg2rad(range(data.shape[0]))
+    [th, r] = np.meshgrid(Theta, rho);
+
+    fig, ax = plt.subplots(figsize=(6,6))
+
+    if ptype == 'mesh':
+        im = ax.pcolormesh(r*np.cos(Theta), r*np.sin(Theta), 100*data.T, 
+                           vmin=vmin,vmax=vmax,cmap=cmap,shading='gouraud')
+    else:
+        im = ax.contourf(r*np.cos(Theta), r*np.sin(Theta), 100*data.T, 
+                           vmin=vmin,vmax=vmax,cmap=cmap,shading='gouraud')
+    cbar = plt.colorbar(im, cax=fig.add_axes([0.15, -0.03, 0.7, 0.05]), orientation='horizontal')
+
+    draw_circle_error(ax)
+    ax.axis('off')
+    if savepath is not None:
+        if not cbar:
+            plt.tight_layout()
+        plt.savefig(savepath, dpi=600)
+    plt.show()
+    
+def plot_bullseye_error(data,mu,vmin=None,vmax=None, savepath=None,cmap='RdBu_r', label='GPRS (%)',n=5):
+    
+    rho     = np.arange(0,4,4.0/data.shape[1]);
+    Theta   = np.deg2rad(range(data.shape[0]))
+    [th, r] = np.meshgrid(Theta, rho);
+
+    fig, ax = plt.subplots(figsize=(6,6))
+
+    levels = np.linspace(vmin, vmax, n+1)
+    im = ax.contourf(r*np.cos(Theta), r*np.sin(Theta), 100*data.T, 
+                           vmin=vmin,vmax=vmax,cmap=cmap,levels=levels)
+
+    cbar = plt.colorbar(im, cax=fig.add_axes([0.15, -0.03, 0.7, 0.05]), orientation='horizontal')
+
+    #ticks = -np.array(range(0,120,20))
+
+    #cbar.set_ticks(ticks);
+    #cbar.set_ticklabels(['%d'%(i) for i in ticks]);
+
+     
+
+    ax.axis('off')
+    draw_circle_error(ax)
+    if savepath is not None:
+        if not cbar:
+            plt.tight_layout()
+        plt.savefig(savepath, dpi=500)
+    plt.show()
+
+def draw_circle_error(ax,width=4):
+
+    circle1 = plt.Circle((0,0), 1, color='black', fill=False, linewidth=width)
+    circle2 = plt.Circle((0,0), 2, color='black', fill=False, linewidth=width)
+    circle3 = plt.Circle((0,0), 3, color='black', fill=False, linewidth=width)
+    circle4 = plt.Circle((0,0), 4, color='black', fill=False, linewidth=width)
+
+    ax.add_artist(circle1)
+    ax.add_artist(circle2)
+    ax.add_artist(circle3)
+    ax.add_artist(circle4)
+    
+    j = 0
+    for i in range(6):
+        theta_i = i*60*np.pi/180 + 60*np.pi/180
+        xi, yi = polar2cart(2, theta_i)
+        xf, yf = polar2cart(4, theta_i)
+        
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+
+    j += 6
+    for i in range(4):
+        theta_i = i*90*(np.pi/180) - 45
+        xi, yi  = polar2cart(1, theta_i)
+        xf, yf  = polar2cart(2, theta_i)
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+
+def draw_circle_frac(ax, mu, width=4, fs=20, color='white'):
+    
+    
+    
+    circle1 = plt.Circle((0,0), 1, color='black', fill=False, linewidth=width)
+    circle2 = plt.Circle((0,0), 2, color='black', fill=False, linewidth=width)
+    circle3 = plt.Circle((0,0), 3, color='black', fill=False, linewidth=width)
+    circle4 = plt.Circle((0,0), 4, color='black', fill=False, linewidth=width)
+
+    ax.add_artist(circle1)
+    ax.add_artist(circle2)
+    ax.add_artist(circle3)
+    ax.add_artist(circle4)
+    
+    j = 0
+    for i in range(6):
+        theta_i = i*60*np.pi/180 + 60*np.pi/180
+        xi, yi = polar2cart(2, theta_i)
+        xf, yf = polar2cart(4, theta_i)
+        
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+        
+        xi, yi = polar2cart(3.5, theta_i + 2*np.pi/12)
+        ax.text(xi-.3, yi, '%.2f' %(mu[j]), weight='bold', fontsize=fs, color=color);
+        xi, yi = polar2cart(2.5, theta_i + 2*np.pi/12)
+        ax.text(xi-.3, yi, '%.2f' %(mu[j+6]), weight='bold', fontsize=fs, color=color); j += 1
+        
+    j += 6
+    LABELS = ['ANT', 'SEPT', 'INF', 'LAT']
+    for i in range(4):
+        theta_i = i*90*np.pi/180 - 45 
+        xi, yi = polar2cart(1, theta_i)
+        xf, yf = polar2cart(2, theta_i)
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+        
+        xi, yi = polar2cart(1.5, theta_i + 2*np.pi/8)
+
+        ax.text(xi-.3, yi, '%.2f' %(mu[j]), weight='bold', fontsize=fs, color=color); j += 1;
+        xi, yi = polar2cart(5, theta_i + 2*np.pi/8)
+
+    ax.text(0-.3, 0-.3, '%.2f' %(mu[j]), weight='bold', fontsize=fs, color=color)
+    
+def draw_circle(ax, mu, width=4, fs=15, xshift=0, yshift=0, color='white'):
+    
+    
+    
+    circle1 = plt.Circle((0,0), 1, color='black', fill=False, linewidth=width)
+    circle2 = plt.Circle((0,0), 2, color='black', fill=False, linewidth=width)
+    circle3 = plt.Circle((0,0), 3, color='black', fill=False, linewidth=width)
+    circle4 = plt.Circle((0,0), 4, color='black', fill=False, linewidth=width)
+
+    ax.add_artist(circle1)
+    ax.add_artist(circle2)
+    ax.add_artist(circle3)
+    ax.add_artist(circle4)
+    
+    j = 0
+    for i in range(6):
+        theta_i = i*60*np.pi/180 + 60*np.pi/180
+        xi, yi = polar2cart(2, theta_i)
+        xf, yf = polar2cart(4, theta_i)
+        
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+        
+        xi, yi = polar2cart(3.5, theta_i + 2*np.pi/12)
+        ax.text(xi-.4-xshift, yi-yshift, '%d' %(mu[j]), weight='bold', fontsize=fs, color=color);
+        xi, yi = polar2cart(2.5, theta_i + 2*np.pi/12)
+        ax.text(xi-.4-xshift, yi-yshift, '%d' %(mu[j+6]), weight='bold', fontsize=fs, color=color); j += 1
+        
+    j += 6
+    LABELS = ['ANT', 'SEPT', 'INF', 'LAT']
+    for i in range(4):
+        theta_i = i*90*np.pi/180  + 45*np.pi/180
+        xi, yi = polar2cart(1, theta_i)
+        xf, yf = polar2cart(2, theta_i)
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+        
+        xi, yi = polar2cart(1.5, theta_i + 2*np.pi/8)
+
+        ax.text(xi-.4-xshift, yi-yshift, '%d' %(mu[j]), weight='bold', fontsize=fs, color=color); j += 1;
+        xi, yi = polar2cart(5, theta_i + 2*np.pi/8)
+
+    ax.text(-.4-xshift, 0-yshift, '%d' %(mu[j]), weight='bold', fontsize=fs, color=color)
+  
+
+def draw_circle_group(ax, mu, std, width=4, fs=14, color='white'):
+    
+        
+    circle1 = plt.Circle((0,0), 1, color='black', fill=False, linewidth=width)
+    circle2 = plt.Circle((0,0), 2, color='black', fill=False, linewidth=width)
+    circle3 = plt.Circle((0,0), 3, color='black', fill=False, linewidth=width)
+    circle4 = plt.Circle((0,0), 4, color='black', fill=False, linewidth=width)
+
+    ax.add_artist(circle1)
+    ax.add_artist(circle2)
+    ax.add_artist(circle3)
+    ax.add_artist(circle4)
+    
+    j = 0
+    for i in range(6):
+        theta_i = i*60*np.pi/180 + 60*np.pi/180
+        xi, yi = polar2cart(2, theta_i)
+        xf, yf = polar2cart(4, theta_i)
+        
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+        
+        xi, yi = polar2cart(3.5, theta_i + 2*np.pi/12)
+        ax.text(xi-.6, yi, '%d(%d)' %(mu[j],std[j]), weight='bold', fontsize=fs, color=color);
+        xi, yi = polar2cart(2.5, theta_i + 2*np.pi/12)
+        ax.text(xi-.6, yi, '%d(%d)' %(mu[j+6],std[j+6]), weight='bold', fontsize=fs, color=color); j += 1
+        
+    j += 6
+    LABELS = ['ANT', 'SEPT', 'INF', 'LAT']
+    for i in range(4):
+        theta_i = i*90*np.pi/180 
+        xi, yi = polar2cart(1, theta_i)
+        xf, yf = polar2cart(2, theta_i)
+        l = Line2D([xi,xf], [yi,yf], color='black', linewidth=width)
+        ax.add_line(l)
+        
+        xi, yi = polar2cart(1.5, theta_i + 2*np.pi/8)
+
+        ax.text(xi-.6, yi-0.1, '%d(%d)' %(mu[j],std[j]), weight='bold', fontsize=fs, color=color); j += 1;
+        xi, yi = polar2cart(5, theta_i + 2*np.pi/8)
+
+    ax.text(0-.3, 0-.2, '%d(%d)' %(mu[j],std[j]), weight='bold', fontsize=fs, color=color)
+    
+def polar2cart(r, theta):
+    x = r * np.cos(theta)
+    y = r * np.sin(theta)
+    return x, y
\ No newline at end of file