Python IO.readPoints方法代码示例

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [as 别名]
def writePoints(filename, points, indices = True):
    """Write points as elastix/transformix point file
    
    Arguments:
        filename (str): file name of the elastix point file.
        points (array or str): source of the points.
        indices (bool): write as pixel indices or physical coordiantes
    
    Returns:
        str : file name of the elastix point file
    """

    points = io.readPoints(points);
    #points = points[:,[1,0,2]]; # points in ClearMap (y,x,z) -> permute to (x,y,z)

  
    with open(filename, 'w') as pointfile:
        if indices:
            pointfile.write('index\n')
        else:
            pointfile.write('point\n')
    
        pointfile.write(str(points.shape[0]) + '\n');
        numpy.savetxt(pointfile, points, delimiter = ' ', newline = '\n', fmt = '%.5e')
        pointfile.close();
    
    return filename;

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [as 别名]
def thresholdPoints(points, intensities, threshold = 0, row = 0):
    """Threshold points by intensities"""
    
    points, intensities = io.readPoints((points, intensities));   
            
    if not isinstance(threshold, tuple):
        threshold = (threshold, all);    
    
    if not isinstance(row, tuple):
        row = (row, row);
    
    
    if intensities.ndim > 1:
        i = intensities[:,row[0]];
    else:
        i = intensities;    
    
    iids = numpy.ones(i.shape, dtype = 'bool');
    if not threshold[0] is all:
        iids = numpy.logical_and(iids, i >= threshold[0]);
        
    if intensities.ndim > 1:
        i = intensities[:,row[1]];
    
    if not threshold[1] is all:
        iids = numpy.logical_and(iids, i <= threshold[1]);
    
    return (points[iids, ...], intensities[iids, ...]);

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [as 别名]
def readPointsGroup(filenames, **args):
    """Turn a list of filenames for points into a numpy stack"""
    
    #check if stack already:
    if isinstance(filenames, numpy.ndarray):
        return filenames;
    
    #read the individual files
    group = [];
    for f in filenames:
        data = io.readPoints(f, **args);
        #data = numpy.reshape(data, (1,) + data.shape);
        group.append(data);
    
    return group

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [as 别名]
def voxelize(points, dataSize = None, sink = None, voxelizeParameter = None,  method = 'Spherical', size = (5,5,5), weights = None):
    """Converts a list of points into an volumetric image array
    
    Arguments:
        points (array): point data array
        dataSize (tuple): size of final image
        sink (str, array or None): the location to write or return the resulting voxelization image, if None return array
        voxelizeParameter (dict):
            ========== ==================== ===========================================================
            Name       Type                 Descritption
            ========== ==================== ===========================================================
            *method*   (str or None)        method for voxelization: 'Spherical', 'Rectangular' or 'Pixel'
            *size*     (tuple)              size parameter for the voxelization
            *weights*  (array or None)      weights for each point, None is uniform weights                          
            ========== ==================== ===========================================================      
    Returns:
        (array): volumetric data of smeared out points
    """
    
    if dataSize is None:
        dataSize = tuple(int(math.ceil(points[:,i].max())) for i in range(points.shape[1]));
    elif isinstance(dataSize, basestring):
        dataSize = io.dataSize(dataSize);
    
    points = io.readPoints(points);
        
    if method.lower() == 'spherical':
        if weights is None:
            data = vox.voxelizeSphere(points.astype('float'), dataSize[0], dataSize[1], dataSize[2], size[0], size[1], size[2]);
        else:
            data = vox.voxelizeSphereWithWeights(points.astype('float'), dataSize[0], dataSize[1], dataSize[2], size[0], size[1], size[2], weights);
           
    elif method.lower() == 'rectangular':
        if weights is None:
            data = vox.voxelizeRectangle(points.astype('float'), dataSize[0], dataSize[1], dataSize[2], size[0], size[1], size[2]);
        else:
            data = vox.voxelizeRectangleWithWeights(points.astype('float'), dataSize[0], dataSize[1], dataSize[2], size[0], size[1], size[2], weights);
    
    elif method.lower() == 'pixel':
        data = voxelizePixel(points, dataSize, weights);
        
    else:
        raise RuntimeError('voxelize: mode: %s not supported!' % method);
    
    return io.writeData(sink, data);

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [as 别名]
def flatfieldFromLine(line, xsize):
    """Creates a 2d flat field image from a 1d line of estimated intensities
    
    Arguments:
        lines (array): array of intensities along y axis
        xsize (int): size of image in x dimension
    
    Returns:
        array: full 2d flat field 
    """
    
    line = io.readPoints(line);

    flatfield = numpy.zeros((xsize, line.size));
    for i in range(xsize):
        flatfield[i,:] = line;
    
    return flatfield;

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [as 别名]
def overlayPoints(dataSource, pointSource, sink = None, pointColor = [1,0,0], x = all, y = all, z = all):
    """Overlay points on 3D data and return as color image
    
    Arguments:
        dataSouce (str or array): volumetric image data
        pointSource (str or array): point data to be overlayed on the image data
        pointColor (array): RGB color for the overlayed points
        x, y, z (all or tuple): sub-range specification
    
    Returns:
        (str or array): image overlayed with points
        
    See Also:
        :func:`overlayLabel`
    """
    data = io.readData(dataSource, x = x, y = y, z = z);
    points = io.readPoints(pointSource, x = x, y = y, z = z, shift = True);
    #print data.shape
    
    if not pointColor is None:
        dmax = data.max(); dmin = data.min();
        if dmin == dmax:
            dmax = dmin + 1;
        cimage = numpy.repeat( (data - dmin) / (dmax - dmin), 3);
        cimage = cimage.reshape(data.shape + (3,));    
    
        if data.ndim == 2:
            for p in points: # faster version using voxelize ?
                cimage[p[0], p[1], :] = pointColor;
        elif data.ndim == 3:
            for p in points: # faster version using voxelize ?
                cimage[p[0], p[1], p[2], :] = pointColor;
        else:
            raise RuntimeError('overlayPoints: data dimension %d not suported' % data.ndim);
    
    else:
        cimage = vox.voxelize(points, data.shape, method = 'Pixel');
        cimage = cimage.astype(data.dtype) * data.max();
        data.shape = data.shape + (1,);
        cimage.shape =  cimage.shape + (1,);
        cimage = numpy.concatenate((data, cimage), axis  = 3);
    
    #print cimage.shape    
    return io.writeData(sink, cimage);   

# 需要导入模块: from ClearMap import IO [as 别名]
# 或者: from ClearMap.IO import readPoints [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 readPoints [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 readPoints [as 别名]
def countPointsInRegions(
    points,
    labeledImage=DefaultLabeledImageFile,
    intensities=None,
    intensityRow=0,
    level=None,
    allIds=False,
    sort=True,
    returnIds=True,
    returnCounts=False,
    collapse=None,
):
    global Label

    points = io.readPoints(points)
    intensities = io.readPoints(intensities)
    pointLabels = labelPoints(points, labeledImage, level=level, collapse=collapse)

    if intensities is None:
        ll, cc = numpy.unique(pointLabels, return_counts=True)
        cci = None
    else:
        if intensities.ndim > 1:
            intensities = intensities[:, intensityRow]

        ll, ii, cc = numpy.unique(pointLabels, return_counts=True, return_inverse=True)
        cci = numpy.zeros(ll.shape)
        for i in range(ii.shape[0]):
            cci[ii[i]] += intensities[i]

    if allIds:
        lla = numpy.setdiff1d(Label.ids, ll)
        ll = numpy.hstack((ll, lla))
        cc = numpy.hstack((cc, numpy.zeros(lla.shape, dtype=cc.dtype)))
        if not cci is None:
            cci = numpy.hstack((cci, numpy.zeros(lla.shape, dtype=cc.dtype)))

    # cc = numpy.vstack((ll,cc)).T;
    if sort:
        ii = numpy.argsort(ll)
        cc = cc[ii]
        ll = ll[ii]
        if not cci is None:
            cci = cci[ii]

    if returnIds:
        if cci is None:
            return ll, cc
        else:
            if returnCounts:
                return ll, cc, cci
            else:
                return ll, cci
    else:
        if cci is None:
            return cc
        else:
            if returnCounts:
                return cc, cci
            else:
                return cci