def interpolateBetweenBinaryObjects(obj1, obj2, slices):
    """
    Takes two binary objects and puts slices slices in-between them, each of which
    contains a smooth binary transition between the objects.
    """
    # constants
    temporal_dimension = 3
    # flip second returned binary objects along temporal axis
    slicer = [slice(None) for _ in range(obj1.ndim + 1)]
    slicer[temporal_dimension] = slice(None, None, -1)
    # logical-and combination
    ret = (
        __interpolateBetweenBinaryObjects(obj1, obj2, slices)
        | __interpolateBetweenBinaryObjects(obj2, obj1, slices)[slicer]
    )
    # control step: see if last volume corresponds to obj2
    slicer[temporal_dimension] = slice(-1, None)
    if not scipy.all(scipy.squeeze(ret[slicer]) == obj2.astype(scipy.bool_)):
        raise Exception(
            "Last created object does not correspond to obj2. Difference of {} voxels.".format(
                len(scipy.nonzero(scipy.squeeze(ret[slicer]) & obj2.astype(scipy.bool_))[0])
            )
        )

    return ret

def Au(U,GF,EpsArr,NX,NY,NZ):
    """Returns the result of matrix-vector multiplication
       by the system matrix A=I-GX
    """
    # reshaping input vector into 4-D array
    Uarr=sci.reshape(U,(NX,NY,NZ,3))
    # extended zero-padded arrays
    Uext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Vext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Jext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    JFext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Uext[0:NX,0:NY,0:NZ,:]=Uarr
    # contrast current array
    s=0
    while s<=2:
        Jext[0:NX,0:NY,0:NZ,s]=Uext[0:NX,0:NY,0:NZ,s]*(EpsArr[0:NX,0:NY,0:NZ]-1.0)
        JFext[:,:,:,s]=fft.fftn(sci.squeeze(Jext[:,:,:,s]))
        s=s+1
    Vext[:,:,:,0]=Uext[:,:,:,0]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,0,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,0,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,0,2],JFext[:,:,:,2])))
    Vext[:,:,:,1]=Uext[:,:,:,1]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,1,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,1,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,1,2],JFext[:,:,:,2])))
    Vext[:,:,:,2]=Uext[:,:,:,2]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,2,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,2,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,2,2],JFext[:,:,:,2])))
    # reshaping output into column vector
    V=sci.reshape(Vext[0:NX,0:NY,0:NZ,:],(NX*NY*NZ*3,1))

    return V

 def dB_Phase_Error_vect(self,other):
     """Subtract two bodes dBmag and phase and return a Bode with
     dBmag and phase."""
     outbode=rwkbode()
     outbode.phase=squeeze(self.phase)-squeeze(other.phase)
     outbode.dBmag=self.dBmag()-other.dBmag()
     return outbode

 def __mul__(self,other):
     """Multiply two bodes."""
     if type(other)==float or type(other)==int:
         myoutput=copy.deepcopy(self)
         myoutput.mag=myoutput.mag*other
         return myoutput
     myin='in'
     myout='out'
     match=1
     if self.output==other.input:
         first=self
         second=other
     elif self.input==other.output:
         first=other
         second=self
     else:
         warnme=1
         if (self.input=='in' and self.output=='out') or (other.input=='in' and other.output=='out'):
             warnme=0
         if warnme:
             print('Warning: multiplying Bodes without a matching input/output pair:\n'+self.output+'/'+self.input+ ' * ' +other.output+'/'+other.input)
         match=0
         first=self
         second=other
     if match:
         myin=first.input
         myout=first.output
     myoutput=copy.deepcopy(self)
     myoutput.input=myin
     myoutput.output=myout
     myoutput.mag=squeeze(colwise(first.mag)*colwise(second.mag))
     myoutput.phase=squeeze(colwise(first.phase)+colwise(second.phase))
     return myoutput

def __load_images(label_image_n, mask_image_n, boundary_image_n, fg_image_n = False, bg_image_n = False):
    """
    Load and return all image data in preprocessed ndarrays.
    The label image will be relabeled to start from 1.
    @return label, ground-truth-mask, boundary-result-mask[, fg-markers, bg-markers]
    """
    # load images
    label_image = load(label_image_n)
    mask_image = load(mask_image_n)
    boundary_image = load(boundary_image_n)
    if fg_image_n: fg_image = load(fg_image_n)
    if bg_image_n: bg_image = load(bg_image_n)
    
    # extract image data
    label_image_d = scipy.squeeze(label_image.get_data())
    mask_image_d = scipy.squeeze(mask_image.get_data()).astype(scipy.bool_)
    boundary_image_d = scipy.squeeze(boundary_image.get_data()).astype(scipy.bool_)
    if fg_image_n: fg_image_d = scipy.squeeze(fg_image.get_data()).astype(scipy.bool_)
    if bg_image_n: bg_image_d = scipy.squeeze(bg_image.get_data()).astype(scipy.bool_)
    
    # check if images are of same dimensionality
    if label_image_d.shape != mask_image_d.shape: raise argparse.ArgumentError('The mask image {} must be of the same dimensionality as the label image {}.'.format(mask_image_d.shape, label_image_d.shape))
    if label_image_d.shape != boundary_image_d.shape: raise argparse.ArgumentError('The boundary term image {} must be of the same dimensionality as the label image {}.'.format(boundary_image_d.shape, label_image_d.shape))
    if fg_image_n:
        if label_image_d.shape != fg_image_d.shape: raise argparse.ArgumentError('The foreground markers image {} must be of the same dimensionality as the label image {}.'.format(fg_image_d.shape, label_image_d.shape))
    if bg_image_n:
        if label_image_d.shape != bg_image_d.shape: raise argparse.ArgumentError('The background markers image {} must be of the same dimensionality as the label image {}.'.format(bg_image_d.shape, label_image_d.shape))
    
    # relabel the label image to start from 1
    label_image_d = medpy.filter.relabel(label_image_d, 1)
    
    if fg_image_n:
        return label_image_d, mask_image_d, boundary_image_d, fg_image_d, bg_image_d
    else:
        return label_image_d, mask_image_d, boundary_image_d

def CombinedBodes(bodelist):
    """This function seeks to make one combined rwkbode instance from
    a list of them.  This may be useful in the case of having several
    experimental Bode data sets with various targeted frequency ranges
    that need to be treated as one data set."""
    docoh=True
    for cb in bodelist:
        if cb.coh is None:
            docoh=False
            break
    first=1
#    pdb.set_trace()
    bigcoh=[]
    for cb in bodelist:
        if first:
            first=0
            bigmag=colwise(cb.mag)
            bigphase=colwise(cb.phase)
            if docoh:
                bigcoh=colwise(cb.coh)
        else:
            bigmag=r_[bigmag,colwise(cb.mag)]
            bigphase=r_[bigphase,colwise(cb.phase)]
            if docoh:
                bigcoh=r_[bigcoh,colwise(cb.coh)]
    myoutbode=copy.deepcopy(bodelist[0])
    myoutbode.mag=squeeze(bigmag)
    myoutbode.phase=squeeze(bigphase)
    myoutbode.coh=squeeze(bigcoh)
    myoutbode.freqlim=[]
    return myoutbode

 def __init__(self, output='out', input='in', \
              mag=None, phase=None, coh=None, \
              freqlim=[], maglim=[], phaselim=[], \
              averaged='not specified', \
              seedfreq=-1, seedphase=0,
              labels=[], legloc=-1, compin=[]):
     self.output = output
     self.input = input
     if len(compin) > 0:
        if mag is None:
           self.mag = squeeze(colwise(abs(compin)))
        if phase is None:
           self.phase = squeeze(colwise(arctan2(imag(compin),real(compin))*180.0/pi))
     else:
         self.mag = squeeze(mag)
         self.phase = squeeze(phase)
     self.coh = coh
     self.averaged = averaged
     self.seedfreq = seedfreq
     self.seedphase = seedphase
     self.freqlim = freqlim
     self.maglim = maglim
     self.phaselim = phaselim
     self.labels = labels
     self.legloc = legloc

 def dB_Phase_Error_sum(self,other, phaseweight=0.1):
     """Subtract two bodes dBmag and phase and return the sum of
     the squared error."""
     phaseerror = squeeze(self.phase)-squeeze(other.phase)
     e_phase = (phaseerror**2).sum()
     dBmagerror = self.dBmag()-other.dBmag()
     e_dB = (dBmagerror**2).sum()
     return e_dB + e_phase*phaseweight

 def ToComp(self,ind=None):
     if ind!=None:
         x=self.mag[ind]*cos(pi/180.0*self.phase[ind])
         y=self.mag[ind]*sin(pi/180.0*self.phase[ind])
         return x+y*1.0j
     else:
         x=squeeze(self.mag)*squeeze(cos(pi/180.0*self.phase))
         y=squeeze(self.mag)*squeeze(sin(pi/180.0*self.phase))
     return squeeze(x+y*1.0j)

def AveBodeFromSpectra(specin):
    Gxxave=squeeze(specin.Gxx.mean(1))
    Gxyave=squeeze(specin.Gxy.mean(1))
    Gyyave=squeeze(specin.Gyy.mean(1))
    Have=Gxyave/Gxxave
    cohnum=norm2(Gxyave)
    cohden=Gxxave*Gyyave
    mybode=rwkbode(input=specin.input, output=specin.output, \
                   compin=squeeze(Have), coh=squeeze(cohnum/cohden))
    return mybode

def RenormalizeFactor(excfile,gsfile,channel=None,Nsamp=1,O=None,q=None):
    if not type(excfile)==list:
        excfile=[excfile]
    if not type(gsfile)==list:
        gsfile=[gsfile]
    exat=GetAttr(excfile[0])
    gsat=GetAttr(gsfile[0])
    L=exat['L']
    if q==None:
        q=sc.array([exat['qx'],exat['qy']])
    if 'phasex' in exat.keys():
        shift=sc.array([exat['phasex']/2.0,exat['phasey']/2.0])
    else:
        shift=sc.array([exat['phase_shift_x']/2.0,exat['phase_shift_y']/2.0])
    phi=exat['phi']
    neel=exat['neel']
    qx,qy,Sq=GetSq(gsfile)
    kx,ky=sf.fermisea(L,L,shift)
    qidx=ml.find((qx==q[0])*(qy==q[1]))
    if O==None:
        _,O,_,_=GetEigSys(excfile,Nsamp)
    pk=None
    sqq=None
    if channel==None:
        channel=exat['channel']
    if channel=='trans':
        pk=sc.squeeze(sf.phiktrans(kx,ky,q[0]/L,q[1]/L,[phi,neel]))
        sqq=sc.real(0.5*(Sq[0,1,qidx]+Sq[0,2,qidx]))
    elif channel=='long':
        pkup=sc.squeeze(sf.phiklong(kx,ky,q[0]/L,q[1]/L,1,[phi,neel]))
        pkdo=sc.squeeze(sf.phiklong(kx,ky,q[0]/L,q[1]/L,-1,[phi,neel]))
        if (q[0]/L==0.5 and q[1]/L==0.5) or (q[0]/L==0 and q[1]/L==0):
            pk=sc.zeros(2*sc.shape(pkup)[0]+1,complex)
        else:
            pk=sc.zeros(2*sc.shape(pkup)[0],complex)
        pk[0:2*sc.shape(pkup)[0]:2]=pkup
        pk[1:2*sc.shape(pkdo)[0]:2]=pkdo
        if (qx[0]/L==0.5 and q[1]/L==0.5) or (q[0]/L==0 and q[1]/L==0):
            if neel==0:
                pk[-1]=0
            else:
                pk[-1]=sum(neel/sf.omega(kx,ky,[phi,neel]))
        sqq=Sq[0,0,qidx]
    else:
        raise(InputFileError('In file \''+excfile+'\', channel=\''+str(channel)+'\'. Should be \'trans\' or \'long\''))
    sqe=sc.einsum('i,jik,k->j',sc.conj(pk),O,pk)
    out=sc.zeros(Nsamp)
    for n in range(Nsamp):
        if abs(sqq)<1e-6 or abs(sqe[n])<1e-6:
            warnings.warn('Probably ill-defined renormalization, returns 1 for sample {0} out of {1}'.format(n,Nsamp),UserWarning)
            out[n]=1
        else:
            out[n]=sc.real(sqq/sqe[n])
    return out

def GetEigSys(filename,gsfile=None,Nsamp=1,channel=None,wavefile=None,q=None):
    if type(filename)==str:
        filename=[filename]
    hfile=h5py.File(filename[0],'r')
    attr=GetAttr(filename[0])
    if channel==None:
        channel=attr['channel']
    dpath,args=GetStat(filename,Nsamp)
    dat=sc.array(hfile["/rank-1/data-0"])
    hfile.close()
    N=int(sc.shape(dat)[0]/2)
    L=attr['L']
    shift=None
    if 'phasex' in attr.keys():
        shift=[attr['phasex']/2.0,attr['phasey']/2.0]
    else:
        shift=[attr['phase_shift_x']/2.0,attr['phase_shift_y']/2.0]
    H=sc.zeros([Nsamp,N,N],complex)
    O=sc.zeros([Nsamp,N,N],complex)
    E=sc.zeros([Nsamp,N])
    V=sc.zeros([Nsamp,N,N],complex)
    for sample,b in enumerate(args):
        for d in b:
            hfile=h5py.File(dpath[d][0],'r')
            dat=hfile[dpath[d][1]]
            H[sample,:,:]+=dat[0:N,0:2*N:2]+1j*dat[0:N,1:2*N:2]
            O[sample,:,:]+=dat[N:2*N,0:2*N:2]+1j*dat[N:2*N,1:2*N:2]
            hfile.close()
        H[sample,:,:]=0.5*(H[sample,:,:]+sc.conj(H[sample,:,:].T))/len(b)
        O[sample,:,:]=0.5*(O[sample,:,:]+sc.conj(O[sample,:,:].T))/len(b)
    if channel=='groundstate':
        return H
    fs=None
    refstate=sc.zeros(2*L*L)
    refstate[0::2]=1
    if wavefile==None:
        fs=GetFermiSigns(filename[0],refstate,channel=channel)
    else:
        fs=GetFermiSigns(wavefile,refstate,channel=channel)
    for s in range(sc.shape(H)[0]):
        H[s,:,:]=sc.dot(sc.diag(fs),sc.dot(H[s,:,:],sc.diag(fs)))
        O[s,:,:]=sc.dot(sc.diag(fs),sc.dot(O[s,:,:],sc.diag(fs)))
    ren=sc.ones(Nsamp)
    if gsfile!=None:
        ren=RenormalizeFactor(filename,gsfile,Nsamp=1,channel=channel,O=O,q=q)
    print('{0} pair of (H,O) matrices loaded, now diagonalize'.format(sc.shape(H)[0]))
    H=sc.einsum('ijk,i->ijk',H,ren)
    O=sc.einsum('ijk,i->ijk',O,ren)
    for s in range(sc.shape(H)[0]):
        E[s,:],V[s,:,:]=vln.geneigh(sc.squeeze(H[s,:,:]),sc.squeeze(O[s,:,:]))
    print('diagonalization finished')
    return H,O,E,V

    def return_results(self, pores=None, throats=None, **kwargs):
        r"""
        Send results of simulation out the the appropriate locations.

        This is a basic version of the update that simply sends out the main
        result (quantity). More elaborate updates should be subclassed.
        """
        if pores is None:
            pores = self.Ps
        if throats is None:
            throats = self.Ts

        phase_quantity = self._quantity.replace(self._phase.name + '_', '')
        if phase_quantity not in self._phase.props():
            self._phase[phase_quantity] = sp.nan
        self._phase[phase_quantity][pores] = self[self._quantity][pores]
        conn_arr = self._net.find_connected_pores(self.Ts)
        dx = sp.squeeze(sp.diff(self[self._quantity][conn_arr], n=1, axis=1))
        g = self['throat.conductance']
        rate = sp.absolute(g * dx)
        if 'throat.rate' not in self._phase.props():
            self._phase['throat.rate'] = sp.nan
        self._phase['throat.rate'][throats] = rate[throats]
        logger.debug('Results of ' + self.name +
                     ' algorithm have been added to ' + self._phase.name)

def BodeFromSpectra(specin,calccoh=False):
    H=specin.Gxy/specin.Gxx
    if calccoh:
       cohnum=squeeze(norm2(specin.Gxy))
       cohden=squeeze(specin.Gxx*specin.Gyy)
       coh=cohnum/cohden
       ave=True
    else:
       coh=[]
       ave=False
    mybode=rwkbode(input=specin.input, output=specin.output, \
                   compin=squeeze(H), coh=coh, averaged=ave)
    mybode.Gxx=specin.Gxx
    mybode.Gyy=specin.Gyy
    mybode.Gxy=specin.Gxy
    return mybode

def plotVol(volume,pts=15,**kwargs):
    fluxGrid = scipy.squeeze(volume.getFluxGrid()).T

    datain = scipy.zeros((fluxGrid.shape[0],
                          pts,
                          fluxGrid.shape[1]))

    for idx in range(pts):
        datain[:,idx,:] = fluxGrid

        temp = genCylGrid(plasma.eq.getRGrid(),
                      scipy.linspace(0,2*scipy.pi,pts),
                      plasma.eq.getZGrid())
    verticies = genVertsFromPixel(temp)

    hex_type = tvtk.api.tvtk.Hexahedron().cell_type
    temp = temp.reshape((temp.size/3,3))
    verticies = verticies.reshape((verticies.size/8,8))

    sg = tvtk.api.tvtk.UnstructuredGrid(points=temp)
    
    sg.set_cells(hex_type,verticies)
    sg.point_data.scalars = datain.flatten()
    sg.point_data.scalars.name = 'temp'
    
    psi_0 = plasma.eq.getFluxAxis()
    psi_LCFS = plasma.eq.getFluxLCFS()
    psi_min = fluxGrid.min()
    psi_max = fluxGrid.max()
    v1 = (psi_0 - psi_min)/(psi_max - psi_min)
    v2 = (psi_LCFS - psi_min)/(psi_max - psi_min)
    mlab.pipeline.volume(sg)

def  __load(picklefile, label):
    """
    Load a pickled testbed as well as the original and label image for further processing.
    The label image will be relabeled to start with region id 1.
    @param picklefile the testbed pickle file name
    @param label the label image file name
    @return a tuple containing:
        label: the label image data as ndarray
        bounding_boxes: the bounding boxes around the label image regions (Note that the the bounding box of a region with id rid is accessed using bounding_boxes[rid - 1])
        model_fg_ids: the region ids of all regions to create the foreground model from
        model_bg_ids: the region ids of all regions to create the background model from
        eval_ids: the regions to evaluate the regions term on, represented by their ids
        truth_fg: subset of regions from the eval_ids that are foreground according to the ground-truth
        truth_bg:  subset of regions from the eval_ids that are background according to the ground-truth
    """
    # load and preprocess images
    label_image = load(label)
    
    label_image_d = scipy.squeeze(label_image.get_data())
    
    # relabel the label image to start from 1
    label_image_d = medpy.filter.relabel(label_image_d, 1)
    
    # extracting bounding boxes
    bounding_boxes = find_objects(label_image_d)
    
    # load testbed
    with open(picklefile, 'r') as f:
        model_fg_ids = cPickle.load(f)
        model_bg_ids = cPickle.load(f)
        cPickle.load(f) # eval ids
        truth_fg = cPickle.load(f)
        truth_bg = cPickle.load(f)
            
    return label_image_d, label_image, bounding_boxes, model_fg_ids, model_bg_ids, truth_fg, truth_bg

def coarsen_2d_ti(Dmat,di=2):
    '''
    Takes a 2D trainin image and makes it coarser, by constructing a 3D TI 
    based on all coarsened 2D images
    '''
    from scipy import squeeze

    nx, ny, nz = Dmat.shape
    ndim3 = di * di    
    x = np.arange(nx)
    y = np.arange(ny)
    ix = x[0:(nx - di):di]
    iy = y[0:(ny - di):di]
    nx2 = len(ix)
    ny2 = len(iy)
    TI = np.zeros( (nx2, ny2, nz, ndim3))
    l = -1;
    for j in range(di):
        for k in range(di):
            l = l + 1
            
            TI_small = Dmat[(0 + j)::di, (0 + k)::di, 0]
            TI[:, :, 0,l] = TI_small[0:nx2, 0:ny2]
            
    TI=squeeze(TI)
    return TI

def xd_iterator_pool(arr, view, fun, processes = None):
    """
    Same as xd_iterator, but distributes the execution over 'processes' subprocesses.
    """
    # check parameters
    for dim in view:
        if dim >= arr.ndim or 0 > dim or not type(dim) == int: raise AttributeError('Invalid dimension {} in view. Must be int between 0 and {}.'.format(dim, arr.ndim - 1))
    if len(view) >= arr.ndim:
        raise AttributeError('The view should contain less entries than the array dimensionality.')
    
    # prepare worker pool    
    worker_pool = multiprocessing.Pool(processes)
    
    # create list of iterations
    iterations = [[None] if dim in view else range(arr.shape[dim]) for dim in range(arr.ndim)]
    
    # prepare views on array (subvolumes)
    slicers = [[slice(None) if idx is None else slice(idx, idx + 1) for idx in indices] for indices in itertools.product(*iterations)]
    reference_shape = arr[slicers[0]].shape
    subvolumes = [scipy.squeeze(arr[slicer]) for slicer in slicers]
        
    # excecute subprocesses
    result_subvolumes = worker_pool.map(fun, subvolumes)
    
    for slicer, result_subvolume in itertools.izip(slicers, result_subvolumes):
        arr[slicer] = result_subvolume.reshape(reference_shape)

    def plot_reals(self, nr=25, hardcopy=0, hardcopy_filename='reals', nanval=-997799, filternan=1):
        import matplotlib.pyplot as plt
        import matplotlib.gridspec as gridspec
        import numpy as np
        from scipy import squeeze
        
        
        nx, ny, nz = self.par['simulation_grid_size']
        
        nr = np.min((self.par['n_real'], nr))
        nsp = int(np.ceil(np.sqrt(nr)))

        fig = plt.figure(1)
        sp = gridspec.GridSpec(nsp, nsp, wspace=0.1, hspace=0.1)
        plt.set_cmap('hot')
        for i in range(0, nr):
            ax1 = plt.Subplot(fig, sp[i])
            fig.add_subplot(ax1)
            if filternan==1:
                self.sim[i][self.sim[i]==nanval] = np.nan
                
            D=squeeze(np.transpose(self.sim[i]))
            plt.imshow(D, extent=[self.x[0], self.x[-1], self.y[0], self.y[-1]], interpolation='none')
            plt.title("Real %d" % (i + 1))

        fig.suptitle(self.method + ' - ' + self.parameter_filename, fontsize=16)
        if (hardcopy):
            plt.savefig(hardcopy_filename)

        plt.show(block=False)

def getArrayFromImage(image):
    """
    Returns a scipy array created from the supplied itk image.
    
    @note always use this instead of directly the itk.PyBuffer, as that object transposes
          the image axes.
    
    @param image an instance of itk.Image holding the array's data
    @type itk.Image (instance)
    
    @return an array
    @rtype scipy.ndarray
    """
    #convert
    itk_py_converter = itk.PyBuffer[getImageType(image)]
    arr = itk_py_converter.GetArrayFromImage(image)
    
    ############
    # !BUG: WarpITK is a very critical itk wrapper. Everything returned / created by it
    # seems to exist only inside the scope of the current function (at least for 
    # ImageFileReader's image and PyBuffers's scipy array. The returned objects have
    # therefore to be copied once, to survive outside the current scope.
    ############
    arr = arr.copy()
    
    # The itk_py_converter transposes the image dimensions. This has to be countered.
    return scipy.squeeze(scipy.transpose(arr))