Python IO.writePoints方法代码示例

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import writePoints [as 别名]
def classifyPixel(project, source, sink = None, processingDirectory = None, cleanup = True):
  """Run pixel classification in headless moded using a trained project file
  
  Arguments:
    project (str): ilastik project .ilp file
    source (str or array): image source
    sink (str or array or None): image sink
  
  Returns:
    str or array: classified image sink
  """
  
  #generate source image file if source is array
  if isinstance(source, str):
    inpfile = source;
  else:
    #generate temporary file
    if processingDirectory is None:
        processingDirectory = tempfile.mkdtemp();
    
    inpfile = os.path.join(processingDirectory, 'ilastik.npy')
    io.writePoints(inpfile, source.transpose((2,1,0)));
  
  if isinstance(sink, str):
    outfile = sink;
  else:
    #outdir  = tempfile.mkdtemp();
    #outfile = os.path.join(outdir, 'result\d*.tif');
    outfile = tempfile.mktemp('.h5');
  
  ilinp = fileNameToIlastikInput(inpfile);
  ilout = fileNameToIlastikOuput(outfile);

  args = '--project="' + project + '" ' + ilout + ' ' + ilinp;
  run(args);
  
  #clean up 
  if cleanup and processingDirectory is not None:
    shutil.rmtree(processingDirectory);
  
  if not isinstance(sink, str):
    sink = readResultH5(outfile);
    if cleanup:
      os.remove(outfile);
  
  return sink;

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import writePoints [as 别名]
def parallelProcessStack(source, x = all, y = all, z = all, sink = None,
                         processes = 2, chunkSizeMax = 100, chunkSizeMin = 30, chunkOverlap = 15,
                         chunkOptimization = True, chunkOptimizationSize = all, 
                         function = noProcessing, join = joinPoints, verbose = False, **parameter):
    """Parallel process a image stack
    
    Main routine that distributes image processing on paralllel processes.
       
    Arguments:
        source (str): image source
        x,y,z (tuple or all): range specifications
        sink (str or None): destination for the result
        processes (int): number of parallel processes
        chunkSizeMax (int): maximal size of a sub-stack
        chunkSizeMin (int): minial size of a sub-stack
        chunkOverlap (int): minimal sub-stack overlap
        chunkOptimization (bool): optimize chunck sizes to best fit number of processes
        chunkOptimizationSize (bool or all): if True only decrease the chunk size when optimizing
        function (function): the main image processing script
        join (function): the fuction to join the results from the image processing script
        verbose (bool): print information on sub-stack generation
        
    Returns:
        str or array: results of the image processing
    """     
    
    subStacks = calculateSubStacks(source, x = x, y = y, z = z, 
                                   processes = processes, chunkSizeMax = chunkSizeMax, chunkSizeMin = chunkSizeMin, chunkOverlap = chunkOverlap,
                                   chunkOptimization = chunkOptimization, chunkOptimizationSize = chunkOptimizationSize, verbose = verbose);
                                   
    nSubStacks = len(subStacks);
    if verbose:
        print "Number of SubStacks: %d" % nSubStacks;
                                       
    #for i in range(nSubStacks):
    #    self.printSubStackInfo(subStacks[i]);
    
    argdata = [];
    for i in range(nSubStacks):
        argdata.append((function, parameter, subStacks[i], verbose));    
    #print argdata
    
    # process in parallel
    pool = Pool(processes = processes);    
    results = pool.map(_processSubStack, argdata);
    
    #print '=========== results';
    #print results;
        
    #join the results
    results = join(results, subStacks = subStacks, **parameter);
    
    #write / or return 
    return io.writePoints(sink, results);

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import writePoints [as 别名]
def transformPoints(source, sink = None, transformParameterFile = None, transformDirectory = None, indices = True, resultDirectory = None, tmpFile = None):
    """Transform coordinates math:`x` via elastix estimated transformation to :math:`T(x)`

    Note the transformation is from the fixed image coorindates to the moving image coordiantes.
    
    Arguments:
        source (str): source of the points
        sink (str or None): sink for transformed points
        transformParameterFile (str or None): parameter file for the primary transformation, if None, the file is determined from the transformDirectory.
        transformDirectory (str or None): result directory of elastix alignment, if None the transformParameterFile has to be given.
        indices (bool): if True use points as pixel coordinates otherwise spatial coordinates.
        resultDirectory (str or None): elastic result directory
        tmpFile (str or None): file name for the elastix point file.
        
    Returns:
        array or str: array or file name of transformed points
    """
        
    global TransformixBinary;    
    
    checkElastixInitialized();    
    global ElastixSettings;

    if tmpFile == None:
        tmpFile = os.path.join(tempfile.tempdir, 'elastix_input.txt');

    # write text file
    if isinstance(source, basestring):
        
        #check if we have elastix signature                 
        with open(source) as f:
            line = f.readline();
            f.close();
            
            if line[:5] == 'point' or line[:5] != 'index':
                txtfile = source;
            else:                
                points = io.readPoints(source);
                #points = points[:,[1,0,2]];
                txtfile = tmpFile;
                writePoints(txtfile, points); 
    
    elif isinstance(source, numpy.ndarray):
        txtfile = tmpFile;
        #points = source[:,[1,0,2]];
        writePoints(txtfile, source);
        
    else:
        raise RuntimeError('transformPoints: source not string or array!');
    
    
    if resultDirectory == None:
        outdirname = os.path.join(tempfile.tempdir, 'elastix_output');
    else:
        outdirname = resultDirectory;
        
    if not os.path.exists(outdirname):
        os.makedirs(outdirname);
        
    
    if transformParameterFile == None:
        if transformDirectory == None:
            RuntimeError('neither alignment directory and transformation parameter file specified!'); 
        transformparameterdir = transformDirectory
        transformparameterfile = getTransformParameterFile(transformparameterdir);
    else:
        transformparameterdir = os.path.split(transformParameterFile);
        transformparameterdir  = transformparameterdir[0];
        transformparameterfile = transformParameterFile;
    
    #transform
    #make path in parameterfiles absolute
    setPathTransformParameterFiles(transformparameterdir);
    
    #run transformix   
    cmd = TransformixBinary + ' -def ' + txtfile + ' -out ' + outdirname + ' -tp ' + transformparameterfile;
    res = os.system(cmd);
    
    if res != 0:
        raise RuntimeError('failed executing ' + cmd);
    
    
    #read data / file 
    if sink == []:
        return io.path.join(outdirname, 'outputpoints.txt')
    
    else:
        #read coordinates
        transpoints = parseElastixOutputPoints(os.path.join(outdirname, 'outputpoints.txt'), indices = indices);

        #correct x,y,z to y,x,z
        #transpoints = transpoints[:,[1,0,2]];     
        
        #cleanup
        for f in os.listdir(outdirname):
            os.remove(os.path.join(outdirname, f));
        os.rmdir(outdirname)
        
        return io.writePoints(sink, transpoints);

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import writePoints [as 别名]
def flatfieldLineFromRegression(data, sink = None, method = 'polynomial', reverse = None, verbose = False):
    """Create flat field line fit from a list of positions and intensities
    
    The fit is either to be assumed to be a Gaussian:
    
    .. math:
        I(x) = a \\exp^{- (x- x_0)^2 / (2 \\sigma)) + b"
        
    or follows a order 6 radial polynomial
        
    .. math:
        I(x) = a + b (x- x_0)^2 + c (x- x_0)^4 + d (x- x_0)^6
    
    Arguments:
        data (array): intensity data as vector of intensities or (n,2) dim array of positions d=0 and intensities measurements d=1:-1 
        sink (str or None): destination to write the result of the fit
        method (str): method to fit intensity data, 'Gaussian' or 'Polynomial'
        reverse (bool): reverse the line fit after fitting
        verbose (bool): print and plot information for the fit
        
    Returns:
        array: fitted intensities on points
    """
    
    data = io.readPoints(data);

    # split data
    if len(data.shape) == 1:
        x = numpy.arange(0, data.shape[0]);
        y = data;
    elif len(data.shape) == 2:
        x = data[:,0]
        y = data[:,1:-1];
    else:
        raise RuntimeError('flatfieldLineFromRegression: input data not a line or array of x,i data');
    
    #calculate mean of the intensity measurements
    ym = numpy.mean(y, axis = 1);

    if verbose:
        plt.figure()
        for i in range(1,data.shape[1]):
            plt.plot(x, data[:,i]);
        plt.plot(x, ym, 'k');
    
    
    if method.lower() == 'polynomial':
        ## fit r^6
        mean = sum(ym * x)/sum(ym)

        def f(x,m,a,b,c,d):
            return a + b * (x-m)**2 + c * (x-m)**4 + d * (x-m)**6;
        
        popt, pcov = curve_fit(f, x, ym, p0 = (mean, 1, 1, 1, .1));
        m = popt[0]; a = popt[1]; b = popt[2];
        c = popt[3]; d = popt[4];

        if verbose:
            print "polynomial fit: %f + %f (x- %f)^2 + %f (x- %f)^4 + %f (x- %f)^6" % (a, b, m, c, m, d, m);

        def fopt(x):
            return f(x, m = m, a = a, b = b, c = c, d = d);
        
        flt = map(fopt, range(0, int(x[-1])));
    
    else: 
        ## Gaussian fit
        
        mean = sum(ym * x)/sum(ym)
        sigma = sum(ym * (x-mean)**2)/(sum(ym))
        
        def f(x, a, m, s, b):
            return a * numpy.exp(- (x - m)**2 / 2 / s) + b;
            
        
        popt, pcov = curve_fit(f, x, ym, p0 = (1000, mean, sigma, 400));
        a = popt[0]; m = popt[1]; s = popt[2]; b = popt[3];

        if verbose:
            print "Gaussian fit: %f exp(- (x- %f)^2 / (2 %f)) + %f" % (a, m, s, b);

        def fopt(x):
            return f(x, a = a, m = m, s = s, b = b);
    
    if reverse:
        flt.reverse();
    
    if verbose:
        plt.plot(x, flt);
        plt.title('flatfieldLineFromRegression')
    
    return io.writePoints(sink, flt);

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import writePoints [as 别名]
def sequentiallyProcessStack(source, x = all, y = all, z = all, sink = None,
                             chunkSizeMax = 100, chunkSizeMin = 30, chunkOverlap = 15,
                             function = noProcessing, join = joinPoints, verbose = False, **parameter):
    """Sequential image processing on a stack
    
    Main routine that sequentially processes a large image on sub-stacks.
       
    Arguments:
        source (str): image source
        x,y,z (tuple or all): range specifications
        sink (str or None): destination for the result
        processes (int): number of parallel processes
        chunkSizeMax (int): maximal size of a sub-stack
        chunkSizeMin (int): minial size of a sub-stack
        chunkOverlap (int): minimal sub-stack overlap
        chunkOptimization (bool): optimize chunck sizes to best fit number of processes
        chunkOptimizationSize (bool or all): if True only decrease the chunk size when optimizing
        function (function): the main image processing script
        join (function): the fuction to join the results from the image processing script
        verbose (bool): print information on sub-stack generation
        
    Returns:
        str or array: results of the image processing
    """     
    #determine z ranges  
    
    subStacks = calculateSubStacks(source, x = x, y = y, z = z, 
                                   processes = 1, chunkSizeMax = chunkSizeMax, chunkSizeMin = chunkSizeMin, chunkOverlap = chunkOverlap,  
                                   chunkOptimization = False, verbose = verbose);
    
    nSubStacks = len(subStacks);
    #print nSubStacks;    
    
    argdata = [];
    for i in range(nSubStacks):
        argdata.append((function, parameter, subStacks[i], verbose));    
    
    #run sequentially
    results = [];
    for i in range(nSubStacks):
        results.append(_processSubStack(argdata[i]));
    
    #join the results
    results = join(results, subStacks = subStacks, **parameter);
    
    #write / or return 
    return io.writePoints(sink, results);










### Pickle does not like classes:

## sub stack information
#class SubStack(object):
#    """Class containing all info of a sub stack usefull for the image processing and result joining functions"""
#    
#    # sub stack id
#    stackId = None;
#    
#    # number of stacks
#    nStacks = None;
#    
#    # tuple of x,y,z range of this sub stack
#    z = all; 
#    x = all;
#    y = all;
#    
#    #original source
#    source = None;    
#    
#    # tuple of center point of the overlaping regions
#    zCenters = None;
#    
#    # tuple of z indices that would generate full image without overlaps
#    zCenterIndices = None;
#    
#    # tuple of z indices in the sub image as returned by readData that would generate full image without overlaps
#    zSubCenterIndices = None;
#    
#    
#    def __init__(slf, stackId = 0, nStacks = 1, source = None, x = all, y = all, z = all, zCenters = all, zCenterIndices = all):
#        slf.stackId = stackId;
#        slf.nStacks = nStacks;
#        slf.source = source; 
#        slf.x = x;
#        slf.y = y;
#        slf.z = z;
#        slf.zCenters = zCenters;
#        slf.zCenterIndices = zCenterIndices;
#        if not zCenterIndices is all and not z is all:
#            slf.zSubCenterIndices = (c - z[0] for c in zCenterIndices);
#        else:
#            slf.zSubCenterIndices = all;
#.........这里部分代码省略.........