First Commit

IP.py

@@ -0,0 +1,243 @@
+from numpy import linspace, zeros, cumsum, ones, where,meshgrid, shape
+from numpy.random import rand, randn
+from skimage.color import rgb2gray
+from skimage.io import imshow, imread
+from skimage.exposure import equalize_hist as histeq
+from skimage.exposure import cumulative_distribution as cdf
+from matplotlib.colors import ListedColormap as colormap
+#from scipy.misc import imresize
+from skimage.transform import resize as imresize
+import matplotlib.pyplot as plt
+from  matplotlib.pyplot import  xlabel, ylabel,title
+from scipy.ndimage.filters import median_filter as medfil2d
+from skimage.filters.rank import maximum
+from skimage.filters.rank import gradient as grad
+from skimage.filters import threshold_otsu as greythres
+import numpy as np
+import skimage.data as Images
+from scipy.signal import convolve2d as filter2D
+from numpy.fft import fft, fft2, ifft, ifft2, fftshift
+from skimage.transform import rotate, radon, iradon, hough_line_peaks,hough_line
+from skimage.morphology import disk,square,rectangle
+from skimage.restoration import wiener
+from skimage.restoration import deconvolution
+from skimage.morphology import binary_erosion as imerode
+from skimage.morphology import binary_dilation as imdilate
+from skimage.morphology import binary_opening as imopen
+from skimage.morphology import binary_closing  as imclose
+from skimage.morphology import erosion as gs_imerode
+from skimage.morphology import dilation as gs_imdilate
+from skimage.morphology import opening as gs_imopen
+from skimage.morphology import closing  as gs_imclose
+from skimage.morphology import black_tophat  as tophat
+from skimage.transform import warp, AffineTransform,ProjectiveTransform
+from skimage import data, color
+from scipy.interpolate import griddata
+from skimage.filters import roberts, sobel, scharr, prewitt
+from skimage.feature import canny, peak_local_max
+from skimage.transform import hough_ellipse
+from skimage.draw import ellipse_perimeter
+from scipy.io import loadmat
+from pylab import imshow, plot, subplot, subplots, show
+
+def bitget(I,Ival):
+    return (I>>Ival) % 2
+
+def hough(img_bin, theta_res=1, rho_res=1):
+    """
+     Computes the Hough transform of an image
+    """
+    nR,nC = img_bin.shape
+    theta = np.linspace(-90.0, 0.0, np.ceil(90.0/theta_res) + 1.0)
+    theta = np.concatenate((theta, -theta[len(theta)-2::-1]))
+
+    D = np.sqrt((nR - 1)**2 + (nC - 1)**2)
+    q = np.ceil(D/rho_res)
+    nrho = 2*q + 1
+    rho = np.linspace(-q*rho_res, q*rho_res, nrho)
+    H = np.zeros((len(rho), len(theta)))
+    for rowIdx in range(nR):
+        for colIdx in range(nC):
+            if img_bin[rowIdx, colIdx]:
+                for thIdx in range(len(theta)):
+                    rhoVal = colIdx*np.cos(theta[thIdx]*np.pi/180.0) + \
+                             rowIdx*np.sin(theta[thIdx]*np.pi/180)
+                    rhoIdx = np.nonzero(np.abs(rho-rhoVal) == np.min(np.abs(rho-rhoVal)))[0]
+                    H[rhoIdx[0], thIdx] += 1
+    return H, theta,rho
+
+
+
+def myImshow(I0,cmp="gray"):
+    imshow(I0,aspect="auto",cmap=cmp)
+def im2bw(Ig,level):
+    S=np.copy(Ig)
+    S[Ig > level] = 1
+    S[Ig <= level] = 0
+    return(S)
+def dtfuv(m,n):
+    """
+    Computes the frequency matrices, used to construct frequency domain filters
+
+    """
+    m = m  # double size to detal with wrap round effects
+    n = n
+    u=linspace(0,m-1,m)
+    v=linspace(0,n-1,n)
+
+    idx = where(u > m/2)
+    u[idx] = u[idx]-m
+
+    idy = where(v > n/2)
+    v[idy] = v[idy]-n
+
+    V,U = meshgrid(v,u)
+    return (V,U)
+
+def fftfilt(f,H):
+    """
+    This function performs fourier domain filtering. I also to zero padding to prevent wrap around effects
+    output filtered imaged
+    f = input image
+    H =  filter coeffienets
+    """
+    f = f.astype("double")
+    I = fft2(f,(shape(H)[0],shape(H)[1]))
+    I2 = ifft2(H*I)
+    return I2[:shape(f)[0],:shape(f)[1]].real
+
+def phantom (n = 256, p_type = 'Modified Shepp-Logan', ellipses = None):
+        """
+         phantom (n = 256, p_type = 'Modified Shepp-Logan', ellipses = None)
+
+        Create a Shepp-Logan or modified Shepp-Logan phantom.
+
+        A phantom is a known object (either real or purely mathematical)
+        that is used for testing image reconstruction algorithms.  The
+        Shepp-Logan phantom is a popular mathematical model of a cranial
+        slice, made up of a set of ellipses.  This allows rigorous
+        testing of computed tomography (CT) algorithms as it can be
+        analytically transformed with the radon transform (see the
+        function `radon').
+
+        Inputs
+        ------
+        n : The edge length of the square image to be produced.
+
+        p_type : The type of phantom to produce. Either
+          "Modified Shepp-Logan" or "Shepp-Logan".  This is overridden
+          if `ellipses' is also specified.
+
+        ellipses : Custom set of ellipses to use.  These should be in
+          the form
+                [[I, a, b, x0, y0, phi],
+                 [I, a, b, x0, y0, phi],
+                 ...]
+          where each row defines an ellipse.
+          I : Additive intensity of the ellipse.
+          a : Length of the major axis.
+          b : Length of the minor axis.
+          x0 : Horizontal offset of the centre of the ellipse.
+          y0 : Vertical offset of the centre of the ellipse.
+          phi : Counterclockwise rotation of the ellipse in degrees,
+                measured as the angle between the horizontal axis and
+                the ellipse major axis.
+          The image bounding box in the algorithm is [-1, -1], [1, 1],
+          so the values of a, b, x0, y0 should all be specified with
+          respect to this box.
+
+        Output
+        ------
+        P : A phantom image.
+
+        Usage example
+        -------------
+          import matplotlib.pyplot as pl
+          P = phantom ()
+          pl.imshow (P)
+
+        References
+        ----------
+        Shepp, L. A.; Logan, B. F.; Reconstructing Interior Head Tissue
+        from X-Ray Transmissions, IEEE Transactions on Nuclear Science,
+        Feb. 1974, p. 232.
+
+        Toft, P.; "The Radon Transform - Theory and Implementation",
+        Ph.D. thesis, Department of Mathematical Modelling, Technical
+        University of Denmark, June 1996.
+
+        """
+
+        if (ellipses is None):
+                ellipses = _select_phantom (p_type)
+        elif (np.size (ellipses, 1) != 6):
+                raise AssertionError ("Wrong number of columns in user phantom")
+
+        # Blank image
+        p = np.zeros ((n, n))
+
+        # Create the pixel grid
+        ygrid, xgrid = np.mgrid[-1:1:(1j*n), -1:1:(1j*n)]
+
+        for ellip in ellipses:
+                I   = ellip [0]
+                a2  = ellip [1]**2
+                b2  = ellip [2]**2
+                x0  = ellip [3]
+                y0  = ellip [4]
+                phi = ellip [5] * np.pi / 180  # Rotation angle in radians
+
+                # Create the offset x and y values for the grid
+                x = xgrid - x0
+                y = ygrid - y0
+
+                cos_p = np.cos (phi)
+                sin_p = np.sin (phi)
+
+                # Find the pixels within the ellipse
+                locs = (((x * cos_p + y * sin_p)**2) / a2
+              + ((y * cos_p - x * sin_p)**2) / b2) <= 1
+
+                # Add the ellipse intensity to those pixels
+                p [locs] += I
+
+        return p
+
+
+def _select_phantom (name):
+        if (name.lower () == 'shepp-logan'):
+                e = _shepp_logan ()
+        elif (name.lower () == 'modified shepp-logan'):
+                e = _mod_shepp_logan ()
+        else:
+                raise ValueError ("Unknown phantom type: %s" % name)
+
+        return e
+
+
+def _shepp_logan ():
+        #  Standard head phantom, taken from Shepp & Logan
+        return [[   2,   .69,   .92,    0,      0,   0],
+                [-.98, .6624, .8740,    0, -.0184,   0],
+                [-.02, .1100, .3100,  .22,      0, -18],
+                [-.02, .1600, .4100, -.22,      0,  18],
+                [ .01, .2100, .2500,    0,    .35,   0],
+                [ .01, .0460, .0460,    0,     .1,   0],
+                [ .02, .0460, .0460,    0,    -.1,   0],
+                [ .01, .0460, .0230, -.08,  -.605,   0],
+                [ .01, .0230, .0230,    0,  -.606,   0],
+                [ .01, .0230, .0460,  .06,  -.605,   0]]
+
+def _mod_shepp_logan ():
+        #  Modified version of Shepp & Logan's head phantom,
+        #  adjusted to improve contrast.  Taken from Toft.
+        return [[   1,   .69,   .92,    0,      0,   0],
+                [-.80, .6624, .8740,    0, -.0184,   0],
+                [-.20, .1100, .3100,  .22,      0, -18],
+                [-.20, .1600, .4100, -.22,      0,  18],
+                [ .10, .2100, .2500,    0,    .35,   0],
+                [ .10, .0460, .0460,    0,     .1,   0],
+                [ .10, .0460, .0460,    0,    -.1,   0],
+                [ .10, .0460, .0230, -.08,  -.605,   0],
+                [ .10, .0230, .0230,    0,  -.606,   0],
+                [ .10, .0230, .0460,  .06,  -.605,   0]]