[SciPy-dev] Extension to PilUtil for inverse radon

Tisham Dhar tisham@apogee.com...
Mon Feb 18 18:11:43 CST 2008


Hi All,

We have added a small inverse radon transform to Pilutil(in scipy/misc),
it goes along with the forward radon transform already there. The file
is below.

# Functions which need the PIL

import types
import numpy

from numpy import amin, amax, ravel, asarray, cast, arange, \
     ones, newaxis, transpose, mgrid, iscomplexobj, sum, zeros, uint8

import math
import Image
import ImageFilter

__all__ = ['fromimage','toimage','imsave','imread','bytescale',
 
'imrotate','imresize','imshow','imfilter','radon','inv_radon']

# Returns a byte-scaled image
def bytescale(data, cmin=None, cmax=None, high=255, low=0):
    if data.dtype is uint8:
        return data
    high = high - low
    if cmin is None:
        cmin = amin(ravel(data))
    if cmax is None:
        cmax = amax(ravel(data))
    scale = high *1.0 / (cmax-cmin or 1)
    bytedata = ((data*1.0-cmin)*scale + 0.4999).astype(uint8)
    return bytedata + cast[uint8](low)

def imread(name,flatten=0):
    """Read an image file from a filename.

    Optional arguments:

     - flatten (0): if true, the image is flattened by calling
convert('F') on
     the resulting image object.  This flattens the color layers into a
single
     grayscale layer.
    """

    im = Image.open(name)
    return fromimage(im,flatten=flatten)

def imsave(name, arr):
    """Save an array to an image file.
    """
    im = toimage(arr)
    im.save(name)
    return

def fromimage(im, flatten=0):
    """Takes a PIL image and returns a copy of the image in a numpy
container.
    If the image is RGB returns a 3-dimensional array:  arr[:,:,n] is
each channel

    Optional arguments:

    - flatten (0): if true, the image is flattened by calling
convert('F') on
    the image object before extracting the numerical data.  This
flattens the
    color layers into a single grayscale layer.  Note that the supplied
image
    object is NOT modified.
    """
    assert Image.isImageType(im), "Not a PIL image."
    if flatten:
        im = im.convert('F')
    mode = im.mode
    adjust = 0
    if mode == '1':
        im = im.convert(mode='L')
        mode = 'L'
        adjust = 1
    str = im.tostring()
    type = uint8
    if mode == 'F':
        type = numpy.float32
    elif mode == 'I':
        type = numpy.uint32
    elif mode == 'I;16':
        type = numpy.uint16
    arr = numpy.fromstring(str,type)
    shape = list(im.size)
    shape.reverse()
    if mode == 'P':
        arr.shape = shape
        if im.palette.rawmode != 'RGB':
            print "Warning: Image has invalid palette."
            return arr
        pal = numpy.fromstring(im.palette.data,type)
        N = len(pal)
        pal.shape = (int(N/3.0),3)
        return arr, pal
    if mode in ['RGB','YCbCr']:
        shape += [3]
    elif mode in ['CMYK','RGBA']:
        shape += [4]
    arr.shape = shape
    if adjust:
        arr = (arr != 0)
    return arr

_errstr = "Mode is unknown or incompatible with input array shape."
def toimage(arr,high=255,low=0,cmin=None,cmax=None,pal=None,
            mode=None,channel_axis=None):
    """Takes a numpy array and returns a PIL image.  The mode of the
    PIL image depends on the array shape, the pal keyword, and the mode
    keyword.

    For 2-D arrays, if pal is a valid (N,3) byte-array giving the RGB
values
    (from 0 to 255) then mode='P', otherwise mode='L', unless mode is
given
    as 'F' or 'I' in which case a float and/or integer array is made

    For 3-D arrays, the channel_axis argument tells which dimension of
the
      array holds the channel data.
    For 3-D arrays if one of the dimensions is 3, the mode is 'RGB'
      by default or 'YCbCr' if selected.
    if the

    The numpy array must be either 2 dimensional or 3 dimensional.
    """
    data = asarray(arr)
    if iscomplexobj(data):
        raise ValueError, "Cannot convert a complex-valued array."
    shape = list(data.shape)
    valid = len(shape)==2 or ((len(shape)==3) and \
                              ((3 in shape) or (4 in shape)))
    assert valid, "Not a suitable array shape for any mode."
    if len(shape) == 2:
        shape = (shape[1],shape[0]) # columns show up first
        if mode == 'F':
            data32 = data.astype(numpy.float32)
            image = Image.fromstring(mode,shape,data32.tostring())
            return image
        if mode in [None, 'L', 'P']:
            bytedata =
bytescale(data,high=high,low=low,cmin=cmin,cmax=cmax)
            image = Image.fromstring('L',shape,bytedata.tostring())
            if pal is not None:
                image.putpalette(asarray(pal,dtype=uint8).tostring())
                # Becomes a mode='P' automagically.
            elif mode == 'P':  # default gray-scale
                pal = arange(0,256,1,dtype=uint8)[:,newaxis] * \
                      ones((3,),dtype=uint8)[newaxis,:]
                image.putpalette(asarray(pal,dtype=uint8).tostring())
            return image
        if mode == '1':  # high input gives threshold for 1
            bytedata = (data > high)
            image = Image.fromstring('1',shape,bytedata.tostring())
            return image
        if cmin is None:
            cmin = amin(ravel(data))
        if cmax is None:
            cmax = amax(ravel(data))
        data = (data*1.0 - cmin)*(high-low)/(cmax-cmin) + low
        if mode == 'I':
            data32 = data.astype(numpy.uint32)
            image = Image.fromstring(mode,shape,data32.tostring())
        else:
            raise ValueError, _errstr
        return image

    # if here then 3-d array with a 3 or a 4 in the shape length.
    # Check for 3 in datacube shape --- 'RGB' or 'YCbCr'
    if channel_axis is None:
        if (3 in shape):
            ca = numpy.flatnonzero(asarray(shape) == 3)[0]
        else:
            ca = numpy.flatnonzero(asarray(shape) == 4)
            if len(ca):
                ca = ca[0]
            else:
                raise ValueError, "Could not find channel dimension."
    else:
        ca = channel_axis

    numch = shape[ca]
    if numch not in [3,4]:
        raise ValueError, "Channel axis dimension is not valid."

    bytedata = bytescale(data,high=high,low=low,cmin=cmin,cmax=cmax)
    if ca == 2:
        strdata = bytedata.tostring()
        shape = (shape[1],shape[0])
    elif ca == 1:
        strdata = transpose(bytedata,(0,2,1)).tostring()
        shape = (shape[2],shape[0])
    elif ca == 0:
        strdata = transpose(bytedata,(1,2,0)).tostring()
        shape = (shape[2],shape[1])
    if mode is None:
        if numch == 3: mode = 'RGB'
        else: mode = 'RGBA'


    if mode not in ['RGB','RGBA','YCbCr','CMYK']:
        raise ValueError, _errstr

    if mode in ['RGB', 'YCbCr']:
        assert numch == 3, "Invalid array shape for mode."
    if mode in ['RGBA', 'CMYK']:
        assert numch == 4, "Invalid array shape for mode."

    # Here we know data and mode is coorect
    image = Image.fromstring(mode, shape, strdata)
    return image

def imrotate(arr,angle,interp='bilinear'):
    """Rotate an image counter-clockwise by angle degrees.

    Interpolation methods can be:
        'nearest' :  for nearest neighbor
        'bilinear' : for bilinear
        'cubic' or 'bicubic' : for bicubic
    """
    arr = asarray(arr)
    func = {'nearest':0,'bilinear':2,'bicubic':3,'cubic':3}
    im = toimage(arr)
    im = im.rotate(angle,resample=func[interp])
    return fromimage(im)

def imresize(arr,newsize,interp='bilinear',mode=None):
    newsize=list(newsize)
    newsize.reverse()
    newsize = tuple(newsize)
    arr = asarray(arr)
    func = {'nearest':0,'bilinear':2,'bicubic':3,'cubic':3}
    im = toimage(arr,mode=mode)
    im = im.resize(newsize,resample=func[interp])
    return fromimage(im)

def imshow(arr):
    """Simple showing of an image through an external viewer.
    """
    im = toimage(arr)
    if (len(arr.shape) == 3) and (arr.shape[2] == 4):
        try:
            import os
            im.save('/tmp/scipy_imshow.png')
            if os.system("(xv /tmp/scipy_imshow.png; rm -f
/tmp/scipy_imshow.png)&"):
                raise RuntimeError
            return
        except:
            print "Warning: Alpha channel may not be handled correctly."

    im.show()
    return

def imresize(arr,size):
    """Resize an image.

    If size is an integer it is a percentage of current size.
    If size is a float it is a fraction of current size.
    If size is a tuple it is the size of the output image.
    """
    im = toimage(arr)
    ts = type(size)
    if ts is types.IntType:
        size = size / 100.0
    if type(size) is types.FloatType:
        size = (im.size[0]*size,im.size[1]*size)
    else:
        size = (size[1],size[0])
    imnew = im.resize(size)
    return fromimage(imnew)


def imfilter(arr,ftype):
    """Simple filtering of an image.

    type can be:
            'blur', 'contour', 'detail', 'edge_enhance',
'edge_enhance_more',
            'emboss', 'find_edges', 'smooth', 'smooth_more', 'sharpen'
    """
    _tdict = {'blur':ImageFilter.BLUR,
              'contour':ImageFilter.CONTOUR,
              'detail':ImageFilter.DETAIL,
              'edge_enhance':ImageFilter.EDGE_ENHANCE,
              'edge_enhance_more':ImageFilter.EDGE_ENHANCE_MORE,
              'emboss':ImageFilter.EMBOSS,
              'find_edges':ImageFilter.FIND_EDGES,
              'smooth':ImageFilter.SMOOTH,
              'smooth_more':ImageFilter.SMOOTH_MORE,
              'sharpen':ImageFilter.SHARPEN
              }

    im = toimage(arr)
    if ftype not in _tdict.keys():
        raise ValueError, "Unknown filter type."
    return fromimage(im.filter(_tdict[ftype]))


def radon(arr,theta=None):
    if theta is None:
        theta = mgrid[0:180]
    s = zeros((arr.shape[1],len(theta)), float)
    k = 0
    for th in theta:
        im = imrotate(arr,-th)
        s[:,k] = sum(im,axis=0)
        k += 1
    return s

def inv_radon(arr):
    theta = mgrid[0:arr.shape[1]]
    th = (math.pi/180)*theta

    #set up the image
    m = arr.shape[1]
    n = arr.shape[0]-1
    inv = zeros( (n, n), int)
    
    #set up x and y matices
    midindex = (n+1)/2
    xy = mgrid[1:n+1,1:n+1]
    xpr = xy[1]-(n+1)/2
    ypr = xy[0]-(n+1)/2

    for l in theta:
        filtrindex = numpy.array(midindex +
xpr*math.sin(th[l])-ypr*math.cos(th[l])).round()
        inva=zeros((n,n))
        spota=numpy.where((filtrindex > 0) & (filtrindex <= n))
        newfiltindex=numpy.array(filtrindex[spota],int)   
        inva[spota]=arr[newfiltindex[:],l]
        inv=inv+inva
    return inv/m

Regards,

Tishampati Dhar

 Remote Sensing Software Developer

 APOGEE IMAGING INTERNATIONAL

 Building 12B

 1 Adelaide - Lobethal Road 

 Lobethal SA 5241 

  Telephone: +61 - 8 - 8389 5499

 Fax: +61 - 8 - 8389 5488

 Mobile: +61 - 406114165

Email: tisham@apogee.com.au mailto:tisham@apogee.com.au>

Web: www.apogee.com.au <http://www.apogee.com.au> 

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

"The information in this e-mail may be confidential and/or commercially
privileged. It is intended solely for the addressee. Access to this
e-mail by anyone else is unauthorised. If you are not the
intendedrecipient, any disclosure, copying, distribution or action taken
or omitted to be taken in reliance on it, is prohibited and may be
unlawful."



More information about the Scipy-dev mailing list