How to use orientation method in Airtest

1# -*- coding: utf-8 -*-2"""3The *Resampling* module provides methods to resample and reorient volumetric 4and point data. 5Resampling the data is usually necessary as the first step to match the 6resolution and orientation of the reference object. 7Main routines for resampling are: :func:`~ClearMap.Alignment.Resampling.resampleData` 8and :func:`~ClearMap.Alignment.Resampling.resamplePoints`.9Image Representation and Size10-----------------------------11The module assumes that images in arrays are arranged as 12    * [x,y] or 13    * [x,y,z] 14where x,y,z correspond to the x,y,z coordinates as displayed in e.g. ImageJ. 15For example an image of size (512,512) stored in an array ``img`` will have:16    >>> img.shape17    (512,512)18Points are assumed to be given as x,y,z coordinates19Parameters such as *resolution* or *dataSize* are assumed to be given in (x,y)20or (x,y,z) format, e.g.21    >>> dataSize = (512,512)22Orientation23-----------24The *orientation* parameter is a tuple of d numbers from 1 to d that specifies 25the permutation of the axes, a minus sign infront of a numbeer indicates26inversion of that axes. For exmaple 27    >>> orientation=(2,-1) 28indicates that x and y should be exchanged and the new y axes should be reversed.29Generally a re-orientation is composed of first a permutation of the axes and then30inverting the indicated axes.31A *permutation* is an orientation without signs and with numbers from 0 to d-1.32Examples:33    >>> import os34    >>> import ClearMap.IO as io35    >>> from ClearMap.Settings import ClearMapPath  36    >>> from ClearMap.Alignment.Resampling import resampleData37    >>> filename = os.path.join(ClearMapPath,'Test/Data/OME/16-17-27_0_8X-s3-20HF_UltraII_C00_xyz-Table Z\d{4}.ome.tif');38    >>> print io.dataSize(filename)39    (2160, 2560, 21)40    >>> data = resampleData(filename, sink = None, resolutionSource = (1,1,1), orientation = (1,2,3), resolutionSink = (10,10,2));41    >>> print data.shape42    (216, 256, 10)43"""44#:copyright: Copyright 2015 by Christoph Kirst, The Rockefeller University, New York City45#:license: GNU, see LICENSE.txt for details.46import sys47import os48import math49import numpy50import pdb51import multiprocessing  52import tempfile53import pdb54import shutil55import cv256#import matplotlib.pyplot as plt57import ClearMap.IO.IO as io58import ClearMap.IO.FileList as fl59from ClearMap.Utils.ProcessWriter import ProcessWriter;60def fixOrientation(orientation):61    """Convert orientation to standard format number sequence62    63    Arguments:64        orientation (tuple or str): orientation specification65        66    Returns:67        tuple: orientation sequence68    69    See Also:70        `Orientation`_71    """72    73    if orientation is None:74        return None;75        76    #fix named representations77    if orientation == 'Left':78        orientation = (1,2,3);79    if orientation == 'Right':80        orientation = (-1,2,3);    81    82    return orientation;83def inverseOrientation(orientation):84    """Returns the inverse permuation of the permutation orientation taking axis inversions into account.85    86    Arguments:87        orientation (tuple or str): orientation specification88        89    Returns:90        tuple: orientation sequence91        92    See Also:93        `Orientation`_94    """95    96    if orientation is None:97        return None;98    99    n = len(orientation);100    iper = list(orientation);101    102    #permutation is defined as permuting the axes and then axis inversion103    for i in range(n):104        if orientation[i] < 0:105            iper[int(abs(orientation[i])-1)] = -(i + 1);106        else:107            iper[int(abs(orientation[i])-1)] = (i + 1);108    109    return tuple(iper)110def orientationToPermuation(orientation):111    """Extracts the permuation from an orientation.112    113    Arguments:114        orientation (tuple or str): orientation specification115        116    Returns:117        tuple: premutation sequence118        119    See Also:120        `Orientation`_121    """122    orientation = fixOrientation(orientation);123    if orientation is None:124        return (0,1,2);125    else:126        return tuple(int(abs(i))-1 for i in orientation);127def orientResolution(resolution, orientation):128    """Permutes a resolution tuple according to the given orientation.129    130    Arguments:131        resolution (tuple): resolution specification132        orientation (tuple or str): orientation specification133        134    Returns:135        tuple: oriented resolution sequence136        137    See Also:138        `Orientation`_139    """140    if resolution is None:141        return None;142    per = orientationToPermuation(orientation);143    #print orientation, per, resolution144    return tuple(resolution[i] for i in per);145    146def orientResolutionInverse(resolution, orientation):147    """Permutes a resolution tuple according to the inverse of a given orientation.148    149    Arguments:150        resolution (tuple): resolution specification151        orientation (tuple or str): orientation specification152        153    Returns:154        tuple: oriented resolution sequence155        156    See Also:157        `Orientation`_158    """159        160    if resolution is None:161        return None;162    163    per = orientationToPermuation(inverseOrientation(orientation));164    return tuple(resolution[i] for i in per);165 166def orientDataSize(dataSize, orientation):167    """Permutes a data size tuple according to the given orientation.168    169    Arguments:170        dataSize (tuple): resolution specification171        orientation (tuple or str): orientation specification172        173    Returns:174        tuple: oriented dataSize sequence175        176    See Also:177        `Orientation`_178    """179    180    return orientResolution(dataSize, orientation);181 182def orientDataSizeInverse(dataSize, orientation):183    """Permutes a dataSize tuple according to the inverse of a given orientation.184    185    Arguments:186        dataSize (tuple): dataSize specification187        orientation (tuple or str): orientation specification188        189    Returns:190        tuple: oriented dataSize sequence191        192    See Also:193        `Orientation`_194    """195    196    return orientResolutionInverse(dataSize, orientation); 197 198 199def resampleDataSize(dataSizeSource, dataSizeSink = None, resolutionSource = None, resolutionSink = None, orientation = None):200    """Calculate scaling factors and data sizes for resampling.201    202    Arguments:203        dataSizeSource (tuple): data size of the original image204        dataSizeSink (tuple or None): data size of the resmapled image205        resolutionSource (tuple or None): resolution of the source image206        resolutionSink (tuple or None): resolution of the sink image207        orientation (tuple or str): re-orientation specification208        209    Returns:210        tuple: data size of the source211        tuple: data size of the sink212        tuple: resolution of source213        tuple: resolution of sink214    215    See Also:216        `Orientation`_217    """218    orientation = fixOrientation(orientation);    219    220    #determine data sizes if not specified221    if dataSizeSink is None:222        if resolutionSource is None or resolutionSink is None:223            raise RuntimeError('resampleDataSize: data size and resolutions not defined!');224        225        #orient resolution of source to resolution of sink to get sink data size226        resolutionSourceO = orientResolution(resolutionSource, orientation);227        dataSizeSourceO = orientDataSize(dataSizeSource, orientation);228        #calculate scaling factor229        dataSizeSink = tuple([int(math.ceil(dataSizeSourceO[i] *  resolutionSourceO[i]/resolutionSink[i])) for i in range(len(dataSizeSource))]);        230        231    #print dataSizeSink, "ds sink"232    233    if dataSizeSource is None:234        if resolutionSource is None or resolutionSink is None:235            raise RuntimeError('resampleDataSize: data size and resolutions not defined!');236        237        #orient resolution of source to resolution of sink to get sink data size238        resolutionSourceO = orientResolution(resolutionSource, orientation);239        240        #calculate source data size241        dataSizeSource = tuple([int(math.ceil(dataSizeSink[i] *  resolutionSink[i]/resolutionSourceO[i])) for i in range(len(dataSizeSink))]);  242        dataSizeSource = orientDataSizeInverse(dataSizeSource);243        244    #print dataSizeSource, "ds source"245        246    #calculate effecive resolutions247    if resolutionSource is None:248        if resolutionSink is None:249            resolutionSource = (1,1,1);250        else:251            dataSizeSourceO = orientDataSize(dataSizeSource, orientation);252            resolutionSource = tuple(float(dataSizeSink[i]) / dataSizeSourceO[i] * resolutionSink[i] for i in range(len(dataSizeSource)));253            resolutionSource = orientResolutionInverse(resolutionSource, orientation);254    255    #print resolutionSource, "res source sink"256    dataSizeSourceO = orientDataSize(dataSizeSource, orientation);257    258    259    resolutionSourceO = orientResolution(resolutionSource, orientation);260    resolutionSink = tuple(float(dataSizeSourceO[i]) / float(dataSizeSink[i]) * resolutionSourceO[i] for i in range(len(dataSizeSource)));261    262    #print dataSizeSource, dataSizeSink, resolutionSource, resolutionSink 263    264    return dataSizeSource, dataSizeSink, resolutionSource, resolutionSink  265def fixInterpolation(interpolation):266    """Converts interpolation given as string to cv2 interpolation object267    268    Arguments:269        interpolation (str or object): interpolation string or cv2 object270    271    Returns:272        object: cv2 interpolation type273    """274    275    if interpolation == 'nn' or interpolation is None or interpolation == cv2.INTER_NEAREST:276        interpolation = cv2.INTER_NEAREST;277    else:278        interpolation = cv2.INTER_LINEAR;279        280    return interpolation;281        282def resampleXY(source, dataSizeSink, sink = None, interpolation = 'linear', out = sys.stdout, verbose = True):283    """Resample a 2d image slice284    285    This routine is used for resampling a large stack in parallel in xy or xz direction.286    287    Arguments:288        source (str or array): 2d image source289        dataSizeSink (tuple): size of the resmapled image290        sink (str or None): location for the resmapled image291        interpolation (str): interpolation method to use: 'linear' or None (nearest pixel)292        out (stdout): where to write progress information293        vebose (bool): write progress info if true294    295    Returns:296        array or str: resampled data or file name297    """   298    299    #out.write("Input: %s Output: " % (inputFile, soutputFile))300    data = io.readData(source);301    dataSize = data.shape;302    303    #print dataSize, dataSizeSink    304    305    if data.ndim != 2:306        raise RuntimeError('resampleXY: expects 2d image source, found %dd' % data.ndim)307    #print sagittalImageSize;308    309    #dataSizeSink = tuple([int(math.ceil(dataSize[i] *  resolutionSource[i]/resolutionSink[i])) for i in range(2)]);310    if verbose:311        out.write(("resampleData: Imagesize: %d, %d " % (dataSize[0], dataSize[1])) + ("Resampled Imagesize: %d, %d" % (dataSizeSink[0], dataSizeSink[1])))312        #out.write(("resampleData: Imagesize: %d, %d " % dataSize) + ("Resampled Imagesize: %d, %d" % (outputSize[1], outputSize[0])))313    314    # note: cv2.resize reverses x-Y axes315    interpolation = fixInterpolation(interpolation)316    sinkData = cv2.resize(data,  (dataSizeSink[1], dataSizeSink[0]), interpolation = interpolation);317    #sinkData = cv2.resize(data,  outputSize);318    #sinkData = scipy.misc.imresize(sagittalImage, outputImageSize, interp = 'bilinear'); #normalizes images -> not usefull for stacks !319    320    #out.write("resampleData: resized Image size: %d, %d " % sinkData.shape)321    322    return io.writeData(sink, sinkData);323def _resampleXYParallel(arg):324    """Resampling helper function to use for parallel resampling of image slices"""325    326    fileSource = arg[0];327    fileSink = arg[1];328    dataSizeSink = arg[2];329    interpolation = arg[3];330    ii = arg[4];331    nn = arg[5];332    verbose = arg[6];333    334    pw = ProcessWriter(ii);335    if verbose:336        pw.write("resampleData: resampling in XY: image %d / %d" % (ii, nn))337    338    data = numpy.squeeze(io.readData(fileSource, z = ii));339    resampleXY(data, sink = fileSink, dataSizeSink = dataSizeSink, interpolation = interpolation, out = pw, verbose = verbose);340def resampleData(source, sink = None,  orientation = None, dataSizeSink = None, resolutionSource = (4.0625, 4.0625, 3), resolutionSink = (25, 25, 25), 341                 processingDirectory = None, processes = 1, cleanup = True, verbose = True, interpolation = 'linear', **args):342    """Resample data of source in resolution and orientation343    344    Arguments:345        source (str or array): image to be resampled346        sink (str or None): destination of resampled image347        orientation (tuple): orientation specified by permuation and change in sign of (1,2,3)348        dataSizeSink (tuple or None): target size of the resampled image349        resolutionSource (tuple): resolution of the source image (in length per pixel)350        resolutionSink (tuple): resolution of the resampled image (in length per pixel)351        processingDirectory (str or None): directory in which to perform resmapling in parallel, None a temporary directry will be created352        processes (int): number of processes to use for parallel resampling353        cleanup (bool): remove temporary files354        verbose (bool): display progress information355        interpolation (str): method to use for interpolating to the resmapled image356    357    Returns:358        (array or str): data or file name of resampled image359    Notes: 360        * resolutions are assumed to be given for the axes of the intrinsic 361          orientation of the data and reference as when viewed by matplotlib or ImageJ362        * orientation: permuation of 1,2,3 with potential sign, indicating which 363          axes map onto the reference axes, a negative sign indicates reversal 364          of that particular axes365        * only a minimal set of information to detremine the resampling parameter 366          has to be given, e.g. dataSizeSource and dataSizeSink367    """368    orientation = fixOrientation(orientation);369    370    if isinstance(dataSizeSink, basestring):371        dataSizeSink = io.dataSize(dataSizeSink);372    #orient actual resolutions onto reference resolution    373    dataSizeSource = io.dataSize(source);374    dataSizeSource, dataSizeSink, resolutionSource, resolutionSink = resampleDataSize(dataSizeSource = dataSizeSource, dataSizeSink = dataSizeSink, 375                                                                                      resolutionSource = resolutionSource, resolutionSink = resolutionSink, orientation = orientation);376    377    dataSizeSinkI = orientDataSizeInverse(dataSizeSink, orientation);378    379    #print dataSizeSource, dataSizeSink, resolutionSource, resolutionSink, dataSizeSinkI380    381     382    #rescale in x y in parallel383    if processingDirectory == None:384        processingDirectory = tempfile.mkdtemp();     385        386    interpolation = fixInterpolation(interpolation);387     388    nZ = dataSizeSource[2];389    pool = multiprocessing.Pool(processes=processes);390    argdata = [];391    for i in range(nZ):392        argdata.append( (source, os.path.join(processingDirectory, 'resample_%04d.tif' % i), dataSizeSinkI, interpolation, i, nZ, verbose) );  393        #print argdata[i]394    pool.map(_resampleXYParallel, argdata);395    396    #rescale in z397    fn = os.path.join(processingDirectory, 'resample_%04d.tif' % 0);398    data = io.readData(fn);399    zImage = numpy.zeros((dataSizeSinkI[0], dataSizeSinkI[1], nZ), dtype = data.dtype);    400    for i in range(nZ):401        if verbose and i % 10 == 0:402            print "resampleData; reading %d/%d" % (i, nZ);403        fn = os.path.join(processingDirectory, 'resample_%04d.tif' % i);404        zImage[:,:, i] = io.readData(fn);405    406    resampledData = numpy.zeros(dataSizeSinkI, dtype = zImage.dtype);407    for i in range(dataSizeSinkI[0]):408        if verbose and i % 25 == 0:409            print "resampleData: processing %d/%d" % (i, dataSizeSinkI[0])410        #resampledImage[:, iImage ,:] =  scipy.misc.imresize(zImage[:,iImage,:], [resizedZAxisSize, sagittalImageSize[1]] , interp = 'bilinear'); 411        #cv2.resize takes reverse order of sizes !412        resampledData[i ,:, :] =  cv2.resize(zImage[i,:,:], (dataSizeSinkI[2], dataSizeSinkI[1]), interpolation = interpolation);413        #resampledData[i ,:, :] =  cv2.resize(zImage[i,:, :], (dataSize[1], resizedZSize));414    415    #account for using (z,y,x) array representation -> (y,x,z)416    #resampledData = resampledData.transpose([1,2,0]);417    #resampledData = resampledData.transpose([2,1,0]);418    419    if cleanup:420        shutil.rmtree(processingDirectory);421    if not orientation is None:422        423        #reorient424        per = orientationToPermuation(orientation);425        resampledData = resampledData.transpose(per);426    427        #reverse orientation after permuting e.g. (-2,1) brings axis 2 to first axis and we can reorder there428        if orientation[0] < 0:429            resampledData = resampledData[::-1, :, :];430        if orientation[1] < 0:431            resampledData = resampledData[:, ::-1, :]; 432        if orientation[2] < 0:433            resampledData = resampledData[:, :, ::-1];434        435        #bring back from y,x,z to z,y,x436        #resampledImage = resampledImage.transpose([2,0,1]);437    if verbose:438        print "resampleData: resampled data size: " + str(resampledData.shape)  439    440    if sink == []:441        if io.isFileExpression(source):442            sink = os.path.split(source);443            sink = os.path.join(sink[0], 'resample_\d{4}.tif');444        elif isinstance(source, basestring):445            sink = source + '_resample.tif';446        else:447            raise RuntimeError('resampleData: automatic sink naming not supported for non string source!');448    scaleFactor = numpy.true_divide(source.shape,resampledData.shape)449    return io.writeData(sink, resampledData), scaleFactor 450    451    452    453def resampleDataInverse(sink, source = None, dataSizeSource = None, orientation = None, resolutionSource = (4.0625, 4.0625, 3), resolutionSink = (25, 25, 25), 454                        processingDirectory = None, processes = 1, cleanup = True, verbose = True, interpolation = 'linear', **args):455    """Resample data inversely to :func:`resampleData` routine456    457    Arguments:458        sink (str or None): image to be inversly resampled (=sink in :func:`resampleData`)459        source (str or array): destination for inversly resmapled image (=source in :func:`resampleData`)460        dataSizeSource (tuple or None): target size of the resampled image461        orientation (tuple): orientation specified by permuation and change in sign of (1,2,3)462        resolutionSource (tuple): resolution of the source image (in length per pixel)463        resolutionSink (tuple): resolution of the resampled image (in length per pixel)464        processingDirectory (str or None): directory in which to perform resmapling in parallel, None a temporary directry will be created465        processes (int): number of processes to use for parallel resampling466        cleanup (bool): remove temporary files467        verbose (bool): display progress information468        interpolation (str): method to use for interpolating to the resmapled image469    470    Returns:471        (array or str): data or file name of resampled image472    Notes: 473        * resolutions are assumed to be given for the axes of the intrinsic 474          orientation of the data and reference as when viewed by matplotlib or ImageJ475        * orientation: permuation of 1,2,3 with potential sign, indicating which 476          axes map onto the reference axes, a negative sign indicates reversal 477          of that particular axes478        * only a minimal set of information to detremine the resampling parameter 479          has to be given, e.g. dataSizeSource and dataSizeSink480    """    481    482    483    #orientation484    orientation = fixOrientation(orientation);485    486    #assume we can read data fully into memory487    resampledData = io.readData(sink);488    dataSizeSink = resampledData.shape;489    490    if isinstance(dataSizeSource, basestring):491        dataSizeSource = io.dataSize(dataSizeSource);492    dataSizeSource, dataSizeSink, resolutionSource, resolutionSink = resampleDataSize(dataSizeSource = dataSizeSource, dataSizeSink = dataSizeSink, 493                                                                                      resolutionSource = resolutionSource, resolutionSink = resolutionSink, orientation = orientation);494    #print (dataSizeSource, dataSizeSink, resolutionSource, resolutionSink )495    496    dataSizeSinkI = orientDataSizeInverse(dataSizeSink, orientation);497    498    499    #flip axes back and permute inversely500    if not orientation is None:501        if orientation[0] < 0:502            resampledData = resampledData[::-1, :, :];503        if orientation[1] < 0:504            resampledData = resampledData[:, ::-1, :]; 505        if orientation[2] < 0:506            resampledData = resampledData[:, :, ::-1];507        508        #reorient509        peri = inverseOrientation(orientation);510        peri = orientationToPermuation(peri);511        resampledData = resampledData.transpose(peri);512    513    # upscale in z514    interpolation = fixInterpolation(interpolation);515    516    resampledDataXY = numpy.zeros((dataSizeSinkI[0], dataSizeSinkI[1], dataSizeSource[2]), dtype = resampledData.dtype);    517    518    for i in range(dataSizeSinkI[0]):519        if verbose and i % 25 == 0:520            print "resampleDataInverse: processing %d/%d" % (i, dataSizeSinkI[0])521        #cv2.resize takes reverse order of sizes !522        resampledDataXY[i ,:, :] =  cv2.resize(resampledData[i,:,:], (dataSizeSource[2], dataSizeSinkI[1]), interpolation = interpolation);523    # upscale x, y in parallel524    525    if io.isFileExpression(source):526        files = source;527    else:528        if processingDirectory == None:529            processingDirectory = tempfile.mkdtemp();   530        files = os.path.join(sink[0], 'resample_\d{4}.tif');531    532    io.writeData(files, resampledDataXY);533    534    nZ = dataSizeSource[2];535    pool = multiprocessing.Pool(processes=processes);536    argdata = [];537    for i in range(nZ):538        argdata.append( (source, fl.fileExpressionToFileName(files, i), dataSizeSource, interpolation, i, nZ) );  539    pool.map(_resampleXYParallel, argdata);540    541    if io.isFileExpression(source):542        return source;543    else:544        data = io.convertData(files, source);545        546        if cleanup:547            shutil.rmtree(processingDirectory);548        549        return data;550    551def resamplePoints(pointSource, pointSink = None, dataSizeSource = None, dataSizeSink = None, orientation = None, resolutionSource = (4.0625, 4.0625, 3), resolutionSink = (25, 25, 25), **args):552    """Resample Points to map from original data to the coordinates of the resampled image553    554    The resampling of points here corresponds to he resampling of an image in :func:`resampleData`555        556    Arguments:557        pointSource (str or array): image to be resampled558        pointSink (str or None): destination of resampled image559        orientation (tuple): orientation specified by permuation and change in sign of (1,2,3)560        dataSizeSource (str, tuple or None): size of the data source561        dataSizeSink (str, tuple or None): target size of the resampled image562        resolutionSource (tuple): resolution of the source image (in length per pixel)563        resolutionSink (tuple): resolution of the resampled image (in length per pixel)564        565    Returns:566        (array or str): data or file name of resampled points567    Notes: 568        * resolutions are assumed to be given for the axes of the intrinsic 569          orientation of the data and reference as when viewed by matplotlib or ImageJ570        * orientation: permuation of 1,2,3 with potential sign, indicating which 571          axes map onto the reference axes, a negative sign indicates reversal 572          of that particular axes573        * only a minimal set of information to detremine the resampling parameter 574          has to be given, e.g. dataSizeSource and dataSizeSink575    """576    577    #fix (y,x,z) image array representation578    #resolutionSource, resolutionSink = self.fixResolutions(resolutionSource, resolutionSink);579    580    orientation = fixOrientation(orientation);581    #datasize of data source582    if isinstance(dataSizeSource, basestring):583        dataSizeSource = io.dataSize(dataSizeSource);584    #orient actual resolutions onto reference resolution    585    dataSizeSource, dataSizeSink, resolutionSource, resolutionSink = resampleDataSize(dataSizeSource = dataSizeSource, dataSizeSink = dataSizeSink, 586                                                                                      resolutionSource = resolutionSource, resolutionSink = resolutionSink, orientation = orientation);587    points = io.readPoints(pointSource);588    dataSizeSinkI = orientDataSizeInverse(dataSizeSink, orientation);589    #resolutionSinkI = orientResolutionInverse(resolutionSink, orientation);590        591    #scaling factors592    temp = len(resolutionSource)593    scale = [float(dataSizeSource[i]) / float(dataSizeSinkI[i]) for i in range(temp)];594    #print scale595    596    repoints = points.copy();597    for i in range(temp):    598        repoints[:,i] = repoints[:,i] / scale[i];599               600    #permute for non trivial orientation601    if not orientation is None:602        per = orientationToPermuation(orientation);603        repoints = repoints[:,per];604        605        for i in range(temp):606            if orientation[i] < 0:607                repoints[:,i] = dataSizeSink[i] - repoints[:,i];608      609    return io.writePoints(pointSink, repoints);610     611def resamplePointsInverse(pointSource, pointSink = None, dataSizeSource = None, dataSizeSink = None, orientation = None, resolutionSource = (4.0625, 4.0625, 3), resolutionSink = (25, 25, 25), **args):612    """Resample points from the coordinates of the resampled image to the original data613    The resampling of points here corresponds to he resampling of an image in :func:`resampleDataInverse`614        615    Arguments:616        pointSource (str or array): image to be resampled617        pointSink (str or None): destination of resampled image618        orientation (tuple): orientation specified by permuation and change in sign of (1,2,3)619        dataSizeSource (str, tuple or None): size of the data source620        dataSizeSink (str, tuple or None): target size of the resampled image621        resolutionSource (tuple): resolution of the source image (in length per pixel)622        resolutionSink (tuple): resolution of the resampled image (in length per pixel)623        624    Returns:625        (array or str): data or file name of inversely resampled points626    Notes: 627        * resolutions are assumed to be given for the axes of the intrinsic 628          orientation of the data and reference as when viewed by matplotlib or ImageJ629        * orientation: permuation of 1,2,3 with potential sign, indicating which 630          axes map onto the reference axes, a negative sign indicates reversal 631          of that particular axes632        * only a minimal set of information to detremine the resampling parameter 633          has to be given, e.g. dataSizeSource and dataSizeSink634    """635       636    orientation = fixOrientation(orientation);637    638    #datasize of data source639    if isinstance(dataSizeSource, basestring):640        dataSizeSource = io.dataSize(dataSizeSource);641    642    dataSizeSource, dataSizeSink, resolutionSource, resolutionSink = resampleDataSize(dataSizeSource = dataSizeSource, dataSizeSink = dataSizeSink, 643                                                                                      resolutionSource = resolutionSource, resolutionSink = resolutionSink, orientation = orientation);644            645    points = io.readPoints(pointSource);646    647    dataSizeSinkI = orientDataSizeInverse(dataSizeSink, orientation);648    #resolutionSinkI = orientResolutionInverse(resolutionSink, orientation);649        650    #scaling factors651    scale = [float(dataSizeSource[i]) / float(dataSizeSinkI[i]) for i in range(3)];652    #print scale653    rpoints = points.copy();    654    655    #invert axis inversion and permutations    656    if not orientation is None:657        #invert permuation658        iorientation = inverseOrientation(orientation);659        per = orientationToPermuation(iorientation);660        rpoints = rpoints[:,per];661        662        for i in range(3):663            if iorientation[i] < 0:664                rpoints[:,i] = dataSizeSink[i] - rpoints[:,i];665    666    #scale points667    for i in range(3):   668        rpoints[:,i] = rpoints[:,i] * scale[i];    669    670    return io.writePoints(pointSink, rpoints);671def sagittalToCoronalData(source, sink = None):672    """Change from saggital to coronal orientation673     674    Arguments:675        source (str or array): source data to be reoriented676        sink (str or None): destination for reoriented image677    678    Returns:679        str or array: reoriented data680    """681      682    source = io.readData(source);683    d = source.ndim;684    if d < 3:685        raise RuntimeError('sagittalToCoronalData: 3d image required!');686    687    tp = range(d);688    tp[0:3] = [2,0,1];689    source = source.transpose(tp);690    source = source[::-1];691    #source = source[::-1,:,:];692    return io.writeData(sink, source);693def _test():694    """Tests for the Resampling Module"""695    import ClearMap.Alignment.Resampling as self696    reload(self)697    from ClearMap.Settings import ClearMapPath as basedir 698    import iDISCO.IO.IO as io699    import os, numpy700    fn = os.path.join(basedir, 'Test/Data/OME/16-17-27_0_8X-s3-20HF_UltraII_C00_xyz-Table Z\d{4}.ome.tif');701    outfn = os.path.join(basedir, "Test/Data/Resampling/test.mhd")702    703    print "Making resampled stack " + outfn704    print "source datasize %s" % str(io.dataSize(fn));705    data = self.resampleData(fn, sink = None, resolutionSource = (1,1,1), orientation = (1,2,3), resolutionSink = (10,10,2));706    print data.shape707    io.writeData(outfn, data)   708    data = self.resampleData(fn, sink = None, dataSizeSink = (50,70,10), orientation = (1,2,3));709    print data.shape710    io.writeData(outfn, data)   711    dataSizeSource, dataSizeSink, resolutionSource, resolutionSink = self.resampleDataSize(dataSizeSource = (100,200, 303), dataSizeSink = None, 712                                                                                      resolutionSource = (1,1,1), resolutionSink = (5,5,5), orientation = (1,2,3));713    print dataSizeSource, dataSizeSink, resolutionSource, resolutionSink714    715    points = numpy.array([[0,0,0], [1,1,1], io.dataSize(fn)]);716    points = points.astype('float')717    pr = self.resamplePoints(points, dataSizeSource = fn, dataSizeSink = (50,70,10), orientation = (1,2,3))718    print pr719    pri = self.resamplePointsInverse(pr, dataSizeSource = fn, dataSizeSink = (50,70,10), orientation = (-1,2,3))720    print pri721    result = self.resampleDataInverse(outfn, os.path.join(basedir, 'Test/Data/OME/resample_\d{4}.ome.tif'), dataSizeSource = fn);722    print result723if __name__ == "__main__":724    _test();725    726#    727#def dataSize(imageFilePattern):728#    """Determine full size from raw data in (x,y,z) order (not the array format (y,x,z))"""729#730#    if os.path.exists(imageFilePattern): # single file731#        tf = tiff.TiffFile(imageFilePattern);732#        shape = tf.series[0]['shape'];733#        return (shape[1], shape[2], shape[0])734#    735#    imageDirectory, listOfImages = self.readFileList(imageFilePattern);736#    nz = len(listOfImages);737#    738#    if nz == 0:739#        raise RuntimeError("dataSize: no files match: %s" % imageFilePattern);740#    741#    imagefile = os.path.join(imageDirectory, listOfImages[0]);742#    sagittalImage = plt.imread(imagefile);# reads as y-x743#    ...

1# Search.py2#3# A* search for wumpus world navigation in Python.4import Action5import Orientation6class SearchState:7	8	def __init__(self, location, orientation, depth, parent, action):9		self.location = location10		self.orientation = orientation11		self.depth = depth12		self.parent = parent13		self.action = action14		self.heuristic = 015		self.cost = 016		17	def __eq__(self, other):18		if ((self.location == other.location) and (self.orientation == other.orientation)):19			return True20		else:21			return False22class SearchEngine:23	24	def __init__(self):25		self.frontier = []26		self.explored = []27		self.safeLocations = []28		self.nodeCount = 029		30	# These are the main methods:31	# - AddSafeLocation: Tell the search about locations you think are safe; the search only considers safe locations to move through.32	# - RemoveSafeLocation: If you determine a safe location is in fact not safe, then you can remove it from consideration.33	# - FindPath: The main method to call to use search to find a sequence of actions leading from start to goal only through safe locations.34	35	def AddSafeLocation (self, x, y):36		if (not self.SafeLocation(x,y)):37			self.safeLocations.append([x,y])38	39	def RemoveSafeLocation(self, x, y):40		if (self.SafeLocation(x,y)):41			self.safeLocations.remove([x,y])42	# Use search to find sequence of actions from start location/orientation to goal location/orientation.43	# Returns empty action list if not path found (or already at the goal).44	def FindPath (self, startLocation, startOrientation, goalLocation, goalOrientation):45		initialState = SearchState (startLocation, startOrientation, 0, None, Action.CLIMB)46		goalState = SearchState (goalLocation, goalOrientation, 0, None, Action.CLIMB)47		finalState = self.Search (initialState, goalState)48		actionList = []49		# If solution found, retain actions50		if (finalState):51			tmpState = finalState52			while (tmpState.parent):53				actionList.insert(0, tmpState.action)54				tmpState = tmpState.parent55		self.Clear() # deletes entire search tree, including initialState and finalState56		return actionList57	# Main search algorithm. Returns goal state from which you can follow the parent pointers58	# to get the actions in the solution path.59	def Search (self,initialState, goalState):60		self.Clear()61		self.nodeCount = 062		print "Calling search...",63		finalState = self.AStarSearch (initialState, goalState)64		if (finalState):65			print "solution found",66		else:67			print "no solution found",68		print " (" + str(self.nodeCount) + " nodes generated)."69		return finalState70	# Clear the explored and frontier lists71	def Clear (self):72		self.frontier = []73		self.explored = []74	# A* search = uniform cost search using cost = (depth + heuristic)75	def AStarSearch (self, initialState, goalState):76		initialState.heuristic = self.CityBlockDistance(initialState.location, goalState.location)77		initialState.cost = initialState.depth + initialState.heuristic78		self.frontier.append(initialState)79		while (self.frontier):80			state = self.frontier.pop(0)81			if (self.GoalTest (state, goalState)):82				return state83			self.explored.append (state)84			# Try each action: GOFORWARD, TURNLEFT, TURNRIGHT85			for action in [Action.GOFORWARD, Action.TURNLEFT, Action.TURNRIGHT]:86				childState = self.GetChildState (state, action)87				if (childState):88					self.nodeCount += 189					childState.heuristic = self.CityBlockDistance (childState.location, goalState.location)90					childState.cost = childState.depth + childState.heuristic91					if (not self.Visited(childState)):92						self.AddToFrontierInOrder(childState)93					else:94						# Check if childState on frontier, but has lower cost95						# Could do this more efficiently by remembering it from Visited call, but good enough96						for index,tmpState in enumerate(self.frontier):97							if (tmpState == childState):98								if (tmpState.cost > childState.cost):99									# Child state is better, so replace frontier state with child state100									self.frontier[index] = childState101								102		return None # failure103	# True if state location same as goal location, ignoring orientation.104	def GoalTest (self, state, goalState):105		if (state == goalState):106			return True107		else:108			return False109	# Returns new state after applying action to given state. For GOFORWARD, only110	# works if adjacent location exists and is safe. If not, return None.111	def GetChildState (self, state, action):112		childState = None113		if (action == Action.TURNLEFT):114			childState = SearchState (state.location, state.orientation, state.depth + 1, state, Action.TURNLEFT)115			if (state.orientation == Orientation.UP):116				childState.orientation = Orientation.LEFT117			if (state.orientation == Orientation.DOWN):118				childState.orientation = Orientation.RIGHT119			if (state.orientation == Orientation.LEFT):120				childState.orientation = Orientation.DOWN121			if (state.orientation == Orientation.RIGHT):122				childState.orientation = Orientation.UP123		if (action == Action.TURNRIGHT):124			childState = SearchState (state.location, state.orientation, state.depth + 1, state, Action.TURNRIGHT)125			if (state.orientation == Orientation.UP):126				childState.orientation = Orientation.RIGHT127			if (state.orientation == Orientation.DOWN):128				childState.orientation = Orientation.LEFT129			if (state.orientation == Orientation.LEFT):130				childState.orientation = Orientation.UP131			if (state.orientation == Orientation.RIGHT):132				childState.orientation = Orientation.DOWN133		if (action == Action.GOFORWARD):134			x = state.location[0]135			y = state.location[1]136			if (state.orientation == Orientation.UP):137				y += 1138			if (state.orientation == Orientation.DOWN):139				y -= 1140			if (state.orientation == Orientation.LEFT):141				x -= 1142			if (state.orientation == Orientation.RIGHT):143				x += 1144			if (self.SafeLocation(x,y)):145				childState = SearchState ([x,y], state.orientation, state.depth + 1, state, Action.GOFORWARD)146		return childState147	def CityBlockDistance (self, location1, location2):148		return (abs (location1[0] - location2[0]) + abs (location1[1] - location2[1]))149	def SafeLocation (self, x, y):150		if ([x,y] in self.safeLocations):151			return True152		else:153			return False154		155	# Return true if state on explored or frontier lists156	def Visited (self, state):157		if (state in self.explored):158			return True159		if (state in self.frontier):160			return True161		return False162	# Insert state into frontier, keeping it in order by state->cost163	# Among equal-cost states, the new state is put first so that A* performs164	# a DFS (more efficient than BFS) among states with equal costs.165	def AddToFrontierInOrder (self, state):166		inserted = False167		for index,tmpState in enumerate(self.frontier):168			if (tmpState.cost >= state.cost):169				self.frontier.insert(index, state)170				inserted = True171				break172		if (not inserted):...