text:
# SimpleCV Image Object

#load required libraries
from SimpleCV.base import *
from SimpleCV.Detection import Barcode, Corner, HaarFeature, Line, Chessboard
from SimpleCV.Features import FeatureSet
from SimpleCV.Stream import JpegStreamer
from SimpleCV.Font import *
from SimpleCV.Color import *
from SimpleCV.DrawingLayer import *

from numpy import int32
from numpy import uint8
import pygame as pg

class ColorSpace:
    """
    This class is used to encapsulates the color space of a given image.
    This class acts like C/C++ style enumerated type.
    See: http://stackoverflow.com/questions/2122706/detect-color-space-with-opencv
    """
    UNKNOWN = 0
    BGR = 1
    GRAY = 2
    RGB = 3
    HLS = 4
    HSV = 5
    XYZ  = 6

class Image:
    """
    The Image class is the heart of SimpleCV and allows you to convert to and
    from a number of source types with ease.  It also has intelligent buffer
    management, so that modified copies of the Image required for algorithms
    such as edge detection, etc can be cached and reused when appropriate.

    Images are converted into 8-bit, 3-channel images in RGB colorspace.  It will
    automatically handle conversion from other representations into this
    standard format.  If dimensions are passed, an empty image is created.
    """
    width = 0    #width and height in px
    height = 0
    depth = 0
    filename = "" #source filename
    filehandle = "" #filehandle if used
    camera = ""
    _mLayers = []

    _barcodeReader = "" #property for the ZXing barcode reader

    #these are buffer frames for various operations on the image
    _bitmap = ""  #the bitmap (iplimage)  representation of the image
    _matrix = ""  #the matrix (cvmat) representation
    _grayMatrix = "" #the gray scale (cvmat) representation -KAS
    _graybitmap = ""  #a reusable 8-bit grayscale bitmap
    _equalizedgraybitmap = "" #the above bitmap, normalized
    _blobLabel = ""  #the label image for blobbing
    _edgeMap = "" #holding reference for edge map
    _cannyparam = (0, 0) #parameters that created _edgeMap
    _pil = "" #holds a PIL object in buffer
    _numpy = "" #numpy form buffer
    _colorSpace = ColorSpace.UNKNOWN #Colorspace Object
    _pgsurface = ""

    #when we empty the buffers, populate with this:
    _initialized_buffers = {
        "_bitmap": "",
        "_matrix": "",
        "_grayMatrix": "",
        "_graybitmap": "",
        "_equalizedgraybitmap": "",
        "_blobLabel": "",
        "_edgeMap": "",
        "_cannyparam": (0, 0),
        "_pil": "",
        "_numpy": "",
        "_pgsurface": ""}

    #initialize the frame
    #parameters: source designation (filename)
    #todo: handle camera/capture from file cases (detect on file extension)
    def __init__(self, source = None, camera = None, colorSpace = ColorSpace.UNKNOWN):
        """
        The constructor takes a single polymorphic parameter, which it tests
        to see how it should convert into an RGB image.  Supported types include:

        OpenCV: iplImage and cvMat types
        Python Image Library: Image type
        Filename: All opencv supported types (jpg, png, bmp, gif, etc)
        """
        self._mLayers = []
        self.camera = camera
        self._colorSpace = colorSpace

        if (type(source) == tuple):
            source = cv.CreateImage(source, cv.IPL_DEPTH_8U, 3)
            cv.Zero(source)
        if (type(source) == cv.cvmat):
            self._matrix = source
            if((source.step/source.cols)==3): #this is just a guess
                self._colorSpace = ColorSpace.BGR
            elif((souce.step/source.cols)==1):
                self._colorSpace = ColorSpace.BGR
            else:
                self._colorSpace = ColorSpace.UNKNOWN

        elif (type(source) == np.ndarray):  #handle a numpy array conversion
            if (type(source[0, 0]) == np.ndarray): #we have a 3 channel array
                #convert to an iplimage bitmap
                source = source.astype(np.uint8)
                source = source[:, :, ::-1].transpose([1, 0, 2])
                self._bitmap = cv.CreateImageHeader((source.shape[1], source.shape[0]), cv.IPL_DEPTH_8U, 3)
                cv.SetData(self._bitmap, source.tostring(),
                    source.dtype.itemsize * 3 * source.shape[1])
                self._numpy = source
                self._colorSpace = ColorSpace.BGR #this is an educated guess
            else:
                #we have a single channel array, convert to an RGB iplimage

                source = source.astype(np.uint8)
                source = source.transpose([1,0]) #we expect width/height but use col/row
                self._bitmap = cv.CreateImage((source.shape[1], source.shape[0]), cv.IPL_DEPTH_8U, 3)
                channel = cv.CreateImageHeader((source.shape[1], source.shape[0]), cv.IPL_DEPTH_8U, 1)
                #initialize an empty channel bitmap
                cv.SetData(channel, source.tostring(),
                    source.dtype.itemsize * source.shape[1])
                cv.Merge(channel, channel, channel, None, self._bitmap)
                self._colorSpace = ColorSpace.BGR

        elif (type(source) == cv.iplimage):
            if (source.nChannels == 1):
                self._bitmap = cv.CreateImage(cv.GetSize(source), cv.IPL_DEPTH_8U, 3)
                cv.Merge(source, source, source, None, self._bitmap)
                self._colorSpace = ColorSpace.BGR
            else:
                self._bitmap = source
                self._colorSpace = ColorSpace.BGR
        elif (type(source) == type(str()) and source != ''):
            self.filename = source
            self._bitmap = cv.LoadImage(self.filename, iscolor=cv.CV_LOAD_IMAGE_COLOR)
            self._colorSpace = ColorSpace.BGR

        elif (type(source) == pg.Surface):
            self._pgsurface = source
            self._bitmap = cv.CreateImageHeader(self._pgsurface.get_size(), cv.IPL_DEPTH_8U, 3)
            cv.SetData(self._bitmap, pg.image.tostring(self._pgsurface, "RGB"))
            cv.CvtColor(self._bitmap, self._bitmap, cv.CV_RGB2BGR)
            self._colorSpace = ColorSpace.BGR

        elif (PIL_ENABLED and source.__class__.__name__ == "JpegImageFile"):
            self._pil = source
            #from the opencv cookbook
            #http://opencv.willowgarage.com/documentation/python/cookbook.html
            self._bitmap = cv.CreateImageHeader(self._pil.size, cv.IPL_DEPTH_8U, 3)
            cv.SetData(self._bitmap, self._pil.tostring())
            self._colorSpace = ColorSpace.BGR
            cv.CvtColor(self._bitmap, self._bitmap, cv.CV_RGB2BGR)
            #self._bitmap = cv.iplimage(self._bitmap)

        else:
            return None
        #if the caller passes in a colorspace we overide it
        if(colorSpace != ColorSpace.UNKNOWN):
            self._colorSpace = colorSpace

        bm = self.getBitmap()
        self.width = bm.width
        self.height = bm.height
        self.depth = bm.depth

    def getColorSpace(self):
        """
        Returns the value matched in the color space class
        so for instance you would use
        if(image.getColorSpace() == ColorSpace.RGB)

        RETURNS: Integer
        """
        return self._colorSpace

    def isRGB(self):
        """
        Returns Boolean
        """
        return(self._colorSpace==ColorSpace.RGB)

    def isBGR(self):
        """
        Returns Boolean
        """
        return(self._colorSpace==ColorSpace.BGR)

    def isHSV(self):
        """
        Returns Boolean
        """
        return(self._colorSpace==ColorSpace.HSV)

    def isHLS(self):
        """
        Returns Boolean
        """
        return(self._colorSpace==ColorSpace.HLS)

    def isXYZ(self):
        """
        Returns Boolean
        """
        return(self._colorSpace==ColorSpace.XYZ)

    def isGray(self):
        """
        Returns Boolean
        """
        return(self._colorSpace==ColorSpace.GRAY)

    def toRGB(self):
        """
        Converts Image colorspace to RGB

        RETURNS: Image
        """
        retVal = self.getEmpty()
        if( self._colorSpace == ColorSpace.BGR or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_BGR2RGB)
        elif( self._colorSpace == ColorSpace.HSV ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HSV2RGB)
        elif( self._colorSpace == ColorSpace.HLS ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HLS2RGB)
        elif( self._colorSpace == ColorSpace.XYZ ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_XYZ2RGB)
        elif( self._colorSpace == ColorSpace.RGB ):
            retVal = self.getBitmap()
        else:
            warnings.warn("Image.toRGB: There is no supported conversion to RGB colorspace")
            return None
        return Image(retVal, colorSpace=ColorSpace.RGB )

    def toBGR(self):
        """
        Converts image colorspace to BGR

        RETURNS: Image
        """
        retVal = self.getEmpty()
        if( self._colorSpace == ColorSpace.RGB or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_RGB2BGR)
        elif( self._colorSpace == ColorSpace.HSV ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HSV2BGR)
        elif( self._colorSpace == ColorSpace.HLS ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HLS2BGR)
        elif( self._colorSpace == ColorSpace.XYZ ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_XYZ2BGR)
        elif( self._colorSpace == ColorSpace.BGR ):
            retVal = self.getBitmap()
        else:
            warnings.warn("Image.toBGR: There is no supported conversion to BGR colorspace")
            return None
        return Image(retVal, colorSpace = ColorSpace.BGR )

    def toHLS(self):
        """
        Converts image to HLS colorspace

        RETURNS: Image
        """
        retVal = self.getEmpty()
        if( self._colorSpace == ColorSpace.BGR or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_BGR2HLS)
        elif( self._colorSpace == ColorSpace.RGB):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_RGB2HLS)
        elif( self._colorSpace == ColorSpace.HSV ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HSV2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2HLS)
        elif( self._colorSpace == ColorSpace.XYZ ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_XYZ2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2HLS)
        elif( self._colorSpace == ColorSpace.HLS ):
            retVal = self.getBitmap()
        else:
            warnings.warn("Image.toHSL: There is no supported conversion to HSL colorspace")
            return None
        return Image(retVal, colorSpace = ColorSpace.HLS )

    def toHSV(self):
        """
        Converts image to HSV colorspace

        RETURNS: Image
        """
        retVal = self.getEmpty()
        if( self._colorSpace == ColorSpace.BGR or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_BGR2HSV)
        elif( self._colorSpace == ColorSpace.RGB):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_RGB2HSV)
        elif( self._colorSpace == ColorSpace.HLS ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HLS2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2HSV)
        elif( self._colorSpace == ColorSpace.XYZ ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_XYZ2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2HSV)
        elif( self._colorSpace == ColorSpace.HSV ):
            retVal = self.getBitmap()
        else:
            warnings.warn("Image.toHSV: There is no supported conversion to HSV colorspace")
            return None
        return Image(retVal, colorSpace = ColorSpace.HSV )

    def toXYZ(self):
        """
        Converts image to XYZ colorspace

        RETURNS: Image
        """
        retVal = self.getEmpty()
        if( self._colorSpace == ColorSpace.BGR or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_BGR2XYZ)
        elif( self._colorSpace == ColorSpace.RGB):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_RGB2XYZ)
        elif( self._colorSpace == ColorSpace.HLS ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HLS2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2XYZ)
        elif( self._colorSpace == ColorSpace.HSV ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HSV2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2XYZ)
        elif( self._colorSpace == ColorSpace.XYZ ):
            retVal = self.getBitmap()
        else:
            warnings.warn("Image.toXYZ: There is no supported conversion to XYZ colorspace")
            return None
        return Image(retVal, colorSpace=ColorSpace.XYZ )

    def toGray(self):
        """
        Converts image to Grayscale colorspace

        RETURNS: Image
        """
        retVal = self.getEmpty(1)
        if( self._colorSpace == ColorSpace.BGR or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_BGR2GRAY)
        elif( self._colorSpace == ColorSpace.RGB):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_RGB2GRAY)
        elif( self._colorSpace == ColorSpace.HLS ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HLS2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2GRAY)
        elif( self._colorSpace == ColorSpace.HSV ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_HSV2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2GRAY)
        elif( self._colorSpace == ColorSpace.XYZ ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_XYZ2RGB)
            cv.CvtColor(retVal, retVal, cv.CV_RGB2GRAY)
        else:
            warnings.warn("Image.toGray: There is no supported conversion to gray colorspace")
            return None
        return Image(retVal, colorSpace = ColorSpace.GRAY )

    def getEmpty(self, channels = 3):
        """
        Create a new, empty OpenCV bitmap with the specified number of channels (default 3)h
        """

        bitmap = cv.CreateImage(self.size(), cv.IPL_DEPTH_8U, channels)
        cv.SetZero(bitmap)
        return bitmap

    def getBitmap(self):
        """
        Retrieve the bitmap (iplImage) of the Image.  This is useful if you want
        to use functions from OpenCV with SimpleCV's image class
        """
        if (self._bitmap):
            return self._bitmap
        elif (self._matrix):
            self._bitmap = cv.GetImage(self._matrix)

        return self._bitmap

    def getMatrix(self):
        """
        Get the matrix (cvMat) version of the image, required for some OpenCV algorithms
        """
        if (self._matrix):
            return self._matrix
        else:
            self._matrix = cv.GetMat(self.getBitmap()) #convert the bitmap to a matrix
            return self._matrix

    def getPIL(self):
        """
        Get a PIL Image object for use with the Python Image Library
        """
        if (not PIL_ENABLED):
            return None
        if (not self._pil):
            rgbbitmap = self.getEmpty()
            cv.CvtColor(self.getBitmap(), rgbbitmap, cv.CV_BGR2RGB)
            self._pil = pil.fromstring("RGB", self.size(), rgbbitmap.tostring())
        return self._pil

    def getNumpy(self):
        """
        Get a Numpy array of the image in width x height x RGB dimensions
        """
        if self._numpy != "":
            return self._numpy

        self._numpy = np.array(self.getMatrix())[:, :, ::-1].transpose([1, 0, 2])
        return self._numpy

    def _getGrayscaleBitmap(self):
        if (self._graybitmap):
            return self._graybitmap

        self._graybitmap = self.getEmpty(1)
        temp = self.getEmpty(3)
        if( self._colorSpace == ColorSpace.BGR or
                self._colorSpace == ColorSpace.UNKNOWN ):
            cv.CvtColor(self.getBitmap(), self._graybitmap, cv.CV_BGR2GRAY)
        elif( self._colorSpace == ColorSpace.RGB):
            cv.CvtColor(self.getBitmap(), self._graybitmap, cv.CV_RGB2GRAY)
        elif( self._colorSpace == ColorSpace.HLS ):
            cv.CvtColor(self.getBitmap(), temp, cv.CV_HLS2RGB)
            cv.CvtColor(temp, self._graybitmap, cv.CV_RGB2GRAY)
        elif( self._colorSpace == ColorSpace.HSV ):
            cv.CvtColor(self.getBitmap(), temp, cv.CV_HSV2RGB)
            cv.CvtColor(temp, self._graybitmap, cv.CV_RGB2GRAY)
        elif( self._colorSpace == ColorSpace.XYZ ):
            cv.CvtColor(self.getBitmap(), retVal, cv.CV_XYZ2RGB)
            cv.CvtColor(temp, self._graybitmap, cv.CV_RGB2GRAY)
        else:
            warnings.warn("Image._getGrayscaleBitmap: There is no supported conversion to gray colorspace")
            return None
        return self._graybitmap

    def getGrayscaleMatrix(self):
        """
        Returns the intensity grayscale matrix
        """
        if (self._grayMatrix):
            return self._grayMatrix
        else:
            self._grayMatrix = cv.GetMat(self._getGrayscaleBitmap()) #convert the bitmap to a matrix
            return self._grayMatrix

    def _getEqualizedGrayscaleBitmap(self):
        if (self._equalizedgraybitmap):
            return self._equalizedgraybitmap

        self._equalizedgraybitmap = self.getEmpty(1)
        cv.EqualizeHist(self._getGrayscaleBitmap(), self._equalizedgraybitmap)

        return self._equalizedgraybitmap

    def getPGSurface(self):
        """
        Gets the pygame surface.  This is used for rendering the display

        RETURNS: pgsurface
        """
        if (self._pgsurface):
            return self._pgsurface
        else:
            self._pgsurface = pg.image.fromstring(self.getPIL().tostring(), self.size(), "RGB")
            return self._pgsurface

    def save(self, filehandle_or_filename="", mode=""):
        """
        Save the image to the specified filename.  If no filename is provided then
        then it will use the filename the Image was loaded from or the last
        place it was saved to.

        Save will implicitly render the image's layers before saving, but the layers are
        not applied to the Image itself.
        """
        if (not filehandle_or_filename):
            if (self.filename):
                filehandle_or_filename = self.filename
            else:
                filehandle_or_filename = self.filehandle

        if (len(self._mLayers)):
            saveimg = self.applyLayers()
        else:
            saveimg = self

        if (type(filehandle_or_filename) != str):
            fh = filehandle_or_filename

            if (not PIL_ENABLED):
                warnings.warn("You need the python image library to save by filehandle")
                return 0

            if (type(fh) == InstanceType and fh.__class__.__name__ == "JpegStreamer"):
                fh.jpgdata = StringIO()
                saveimg.getPIL().save(fh.jpgdata, "jpeg") #save via PIL to a StringIO handle
                fh.refreshtime = time.time()
                self.filename = ""
                self.filehandle = fh

            elif (type(fh) == InstanceType and fh.__class__.__name__ == "VideoStream"):
                self.filename = ""
                self.filehandle = fh
                fh.writeFrame(saveimg)

            elif (type(fh) == InstanceType and fh.__class__.__name__ == "Display"):
                self.filename = ""
                self.filehandle = fh
                fh.writeFrame(saveimg)

            else:
                if (not mode):
                    mode = "jpeg"

                saveung.getPIL().save(fh, mode)
                self.filehandle = fh #set the filename for future save operations
                self.filename = ""

            return 1

        filename = filehandle_or_filename
        if (filename):
            cv.SaveImage(filename, saveimg.getBitmap())
            self.filename = filename #set the filename for future save operations
            self.filehandle = ""
        elif (self.filename):
            cv.SaveImage(self.filename, saveimg.getBitmap())
        else:
            return 0

        return 1

    def copy(self):
        """
        Return a full copy of the Image's bitmap.  Note that this is different
        from using python's implicit copy function in that only the bitmap itself
        is copied.

        Returns: IMAGE
        """
        newimg = self.getEmpty()
        cv.Copy(self.getBitmap(), newimg)
        return Image(newimg, colorSpace=self._colorSpace)

    #scale this image, and return a new Image object with the new dimensions
    def scale(self, width, height = -1):
        """
        Scale the image to a new width and height.

        If no height is provided, the width is considered a scaling value ie::

            img.scale(200, 100) #scales the image to 200px x 100px
            img.scale(2.0) #enlarges the image to 2x its current size

        Returns: IMAGE
        """
        w, h = width, height
        if height == -1:
          w = int(self.width * width)
          h = int(self.height * width)

        scaled_bitmap = cv.CreateImage((w, h), 8, 3)
        cv.Resize(self.getBitmap(), scaled_bitmap)
        return Image(scaled_bitmap, colorSpace=self._colorSpace)

    def smooth(self, algorithm_name = 'gaussian', aperature = '', sigma = 0, spatial_sigma = 0):
        """
        Smooth the image, by default with the Gaussian blur.  If desired,
        additional algorithms and aperatures can be specified.  Optional parameters
        are passed directly to OpenCV's cv.Smooth() function.

        Returns: IMAGE
        """
        win_x = 3
        win_y = 3  #set the default aperature window size (3x3)

        if (is_tuple(aperature)):
            win_x, win_y = aperature#get the coordinates from parameter
            #TODO: make sure aperature is valid
            #   eg Positive, odd and square for bilateral and median

        algorithm = cv.CV_GAUSSIAN #default algorithm is gaussian

        #gauss and blur can work in-place, others need a buffer frame
        #use a string to ID rather than the openCV constant
        if algorithm_name == "blur":
            algorithm = cv.CV_BLUR
        if algorithm_name == "bilateral":
            algorithm = cv.CV_BILATERAL
            win_y = win_x #aperature must be square
        if algorithm_name == "median":
            algorithm = cv.CV_MEDIAN
            win_y = win_x #aperature must be square

        newimg = self.getEmpty(1)
        cv.Smooth(self._getGrayscaleBitmap(), newimg, algorithm, win_x, win_y, sigma, spatial_sigma)

        return Image(newimg, colorSpace=self._colorSpace)

    def medianFilter(self, window=''):
        """
        Perform a median filtering operation to denoise/despeckle the image.
        The optional parameter is the window size.
        """
        return self.smooth(algorithm_name='median', aperature=window)

    def bilateralFilter(self, window = ''):
        """
        Perform a bilateral filtering operation to denoise/despeckle the image.
        The optional parameter is the window size.
        """
        return self.smooth(algorithm_name='bilateral', aperature=window)

    def invert(self):
        """
        Invert (negative) the image note that this can also be done with the
        unary minus (-) operator.

        Returns: IMAGE
        """
        return -self

    def grayscale(self):
        """
        return a gray scale version of the image

        Returns: IMAGE
        """
        return Image(self._getGrayscaleBitmap())

    def flipHorizontal(self):
        """
        Horizontally mirror an image

        Returns: IMAGE
        """
        newimg = self.getEmpty()
        cv.Flip(self.getBitmap(), newimg, 1)
        return Image(newimg, colorSpace=self._colorSpace)

    def flipVertical(self):
        """
        Vertically mirror an image

        Returns: IMAGE
        """
        newimg = self.getEmpty()
        cv.Flip(self.getBitmap(), newimg, 0)
        return Image(newimg, colorSpace=self._colorSpace)

    def stretch(self, thresh_low = 0, thresh_high = 255):
        """
        The stretch filter works on a greyscale image, if the image
        is color, it returns a greyscale image.  The filter works by
        taking in a lower and upper threshold.  Anything below the lower
        threshold is pushed to black (0) and anything above the upper
        threshold is pushed to white (255)

        Returns: IMAGE
        """
        try:
            newimg = self.getEmpty(1)
            cv.Threshold(self._getGrayscaleBitmap(), newimg, thresh_low, 255, cv.CV_THRESH_TOZERO)
            cv.Not(newimg, newimg)
            cv.Threshold(newimg, newimg, 255 - thresh_high, 255, cv.CV_THRESH_TOZERO)
            cv.Not(newimg, newimg)
            return Image(newimg)
        except:
            return None

    def binarize(self, thresh = 127, maxv = 255, blocksize = 3, p = 5):
        """
        Do a binary threshold the image, changing all values above thresh to maxv
        and all below to black.  If a color tuple is provided, each color channel
        is thresholded separately.

        If threshold is -1, an adaptive sampling method is used - similar to a moving
        average.  Over each region of block*block pixels a threshold is applied
        where threshold = local_mean - p.
        """
        if (is_tuple(thresh)):
            r = self.getEmpty(1)
            g = self.getEmpty(1)
            b = self.getEmpty(1)
            cv.Split(self.getBitmap(), b, g, r, None)

            cv.Threshold(r, r, thresh[0], maxv, cv.CV_THRESH_BINARY)
            cv.Threshold(g, g, thresh[1], maxv, cv.CV_THRESH_BINARY)
            cv.Threshold(b, b, thresh[2], maxv, cv.CV_THRESH_BINARY)

            cv.Add(r, g, r)
            cv.Add(r, b, r)

            return Image(r, colorSpace=self._colorSpace)

        elif thresh == -1:
            newbitmap = self.getEmpty(1)
            cv.AdaptiveThreshold(self._getGrayscaleBitmap(), newbitmap, maxv,
                cv.CV_ADAPTIVE_THRESH_GAUSSIAN_C, cv.CV_THRESH_BINARY_INV, blocksize, p)
            return Image(newbitmap, colorSpace=self._colorSpace)
        else:
            newbitmap = self.getEmpty(1)
            #desaturate the image, and apply the new threshold
            cv.Threshold(self._getGrayscaleBitmap(), newbitmap, thresh, float(maxv), cv.CV_THRESH_BINARY_INV)
            return Image(newbitmap, colorSpace=self._colorSpace)

    def meanColor(self):
        """
        Finds average color of all the pixels in the image.

        Returns: IMAGE
        """
        return tuple(reversed(cv.Avg(self.getBitmap())[0:3]))

    def findCorners(self, maxnum = 50, minquality = 0.04, mindistance = 1.0):
        """
        This will find corner Feature objects and return them as a FeatureSet
        strongest corners first.  The parameters give the number of corners to look
        for, the minimum quality of the corner feature, and the minimum distance
        between corners.

        Returns: FEATURESET
        """
        #initialize buffer frames
        eig_image = cv.CreateImage(cv.GetSize(self.getBitmap()), cv.IPL_DEPTH_32F, 1)
        temp_image = cv.CreateImage(cv.GetSize(self.getBitmap()), cv.IPL_DEPTH_32F, 1)

        corner_coordinates = cv.GoodFeaturesToTrack(self._getGrayscaleBitmap(), eig_image, temp_image, maxnum, minquality, mindistance, None)

        corner_features = []
        for (x, y) in corner_coordinates:
            corner_features.append(Corner(self, x, y))

        return FeatureSet(corner_features)

    def findBlobs(self, threshval = 127, minsize=10, maxsize=0, threshblocksize=3, threshconstant=5):
        """
        This will look for continuous
        light regions and return them as Blob features in a FeatureSet.  Parameters
        specify the binarize filter threshold value, and minimum and maximum size for blobs.
        If a threshold value is -1, it will use an adaptive threshold.  See binarize() for
        more information about thresholding.  The threshblocksize and threshconstant
        parameters are only used for adaptive threshold.

        Note that this previously used cvblob and the python-cvblob library,
        which is no longer necessary

        Returns: FEATURESET
        """
        if (maxsize == 0):
            maxsize = self.width * self.height / 2
        #create a single channel image, thresholded to parameters

        blobmaker = BlobMaker()
        blobs = blobmaker.extractFromBinary(self.binarize(threshval, 255, threshblocksize, threshconstant).invert(),
            self, minsize = minsize, maxsize = maxsize)

        if not len(blobs):
            return None

        return FeatureSet(blobs).sortArea()

    #this code is based on code that's based on code from
    #http://blog.jozilla.net/2008/06/27/fun-with-python-opencv-and-face-detection/
    def findHaarFeatures(self, cascade, scale_factor=1.2, min_neighbors=2, use_canny=cv.CV_HAAR_DO_CANNY_PRUNING):
        """
        If you want to find Haar Features (useful for face detection among other
        purposes) this will return Haar feature objects in a FeatureSet.  The
        parameters are:
        * the scaling factor for subsequent rounds of the haar cascade (default 1.2)7
        * the minimum number of rectangles that makes up an object (default 2)
        * whether or not to use Canny pruning to reject areas with too many edges (default yes, set to 0 to disable)

        For more information, consult the cv.HaarDetectObjects documentation

        You will need to provide your own cascade file - these are usually found in
        /usr/local/share/opencv/haarcascades and specify a number of body parts.

        Note that the cascade parameter can be either a filename, or a HaarCascade
        loaded with cv.Load().

        Returns: FEATURESET
        """
        storage = cv.CreateMemStorage(0)

        #lovely.  This segfaults if not present
        if type(cascade) == str:
          if (not os.path.exists(cascade)):
              warnings.warn("Could not find Haar Cascade file " + cascade)
              return None
          cascade = cv.Load(cascade)

        objects = cv.HaarDetectObjects(self._getEqualizedGrayscaleBitmap(), cascade, storage, scale_factor, use_canny)
        if objects:
            return FeatureSet([HaarFeature(self, o, cascade) for o in objects])

        return None

    def drawCircle(self, ctr, rad, color = (0, 0, 0), thickness = 1):
        """
        Draw a circle on the Image, parameters include:
        * the center of the circle
        * the radius in pixels
        * a color tuple (default black)
        * the thickness of the circle

        Note that this function is depricated, try to use DrawingLayer.circle() instead

        Returns: NONE - Inline Operation
        """
        self.getDrawingLayer().circle((int(ctr[0]), int(ctr[1])), int(rad), color, int(thickness))

    def drawLine(self, pt1, pt2, color = (0, 0, 0), thickness = 1):
        """
        Draw a line on the Image, parameters include
        * pt1 - the first point for the line (tuple)
        * pt1 - the second point on the line (tuple)
        * a color tuple (default black)
        * thickness of the line

        Note that this modifies the image in-place and clears all buffers.

        Returns: NONE - Inline Operation
        """
        pt1 = (int(pt1[0]), int(pt1[1]))
        pt2 = (int(pt2[0]), int(pt2[1]))
        self.getDrawingLayer().line(pt1, pt2, color, thickness)

    def size(self):
        """
        Gets width and height

        Returns: TUPLE
        """
        return cv.GetSize(self.getBitmap())

    def split(self, cols, rows):
        """
        Given number of cols and rows, splits the image into a cols x rows 2d array
        of cropped images

        quadrants = Image("foo.jpg").split(2,2) <-- returns a 2d array of 4 images
        """
        crops = []

        wratio = self.width / cols
        hratio = self.height / rows

        for i in range(rows):
            row = []
            for j in range(cols):
                row.append(self.crop(j * wratio, i * hratio, wratio, hratio))
            crops.append(row)

        return crops

    def splitChannels(self, grayscale = True):
        """
        Split the channels of an image into RGB (not the default BGR)
        single parameter is whether to return the channels as grey images (default)
        or to return them as tinted color image

        Returns: TUPLE - of 3 image objects
        """
        r = self.getEmpty(1)
        g = self.getEmpty(1)
        b = self.getEmpty(1)
        cv.Split(self.getBitmap(), b, g, r, None)

        red = self.getEmpty()
        green = self.getEmpty()
        blue = self.getEmpty()

        if (grayscale):
            cv.Merge(r, r, r, None, red)
            cv.Merge(g, g, g, None, green)
            cv.Merge(b, b, b, None, blue)
        else:
            cv.Merge(None, None, r, None, red)
            cv.Merge(None, g, None, None, green)
            cv.Merge(b, None, None, None, blue)

        return (Image(red), Image(green), Image(blue))

    def applyHLSCurve(self, hCurve, lCurve, sCurve):
        """
        Apply 3 ColorCurve corrections applied in HSL space
        Parameters are:
        * Hue ColorCurve
        * Lightness (brightness/value) ColorCurve
        * Saturation ColorCurve

        Returns: IMAGE
        """

        #TODO CHECK ROI
        #TODO CHECK CURVE SIZE
        #TODO CHECK COLORSPACE
        #TODO CHECK CURVE SIZE
        temp  = cv.CreateImage(self.size(), 8, 3)
        #Move to HLS space
        cv.CvtColor(self._bitmap, temp, cv.CV_RGB2HLS)
        tempMat = cv.GetMat(temp) #convert the bitmap to a matrix
        #now apply the color curve correction
        tempMat = np.array(self.getMatrix()).copy()
        tempMat[:, :, 0] = np.take(hCurve.mCurve, tempMat[:, :, 0])
        tempMat[:, :, 1] = np.take(sCurve.mCurve, tempMat[:, :, 1])
        tempMat[:, :, 2] = np.take(lCurve.mCurve, tempMat[:, :, 2])
        #Now we jimmy the np array into a cvMat
        image = cv.CreateImageHeader((tempMat.shape[1], tempMat.shape[0]), cv.IPL_DEPTH_8U, 3)
        cv.SetData(image, tempMat.tostring(), tempMat.dtype.itemsize * 3 * tempMat.shape[1])
        cv.CvtColor(image, image, cv.CV_HLS2RGB)
        return Image(image, colorSpace=self._colorSpace)

    def applyRGBCurve(self, rCurve, gCurve, bCurve):
        """
        Apply 3 ColorCurve corrections applied in rgb channels
        Parameters are:
        * Red ColorCurve
        * Green ColorCurve
        * Blue ColorCurve

        Returns: IMAGE
        """
        tempMat = np.array(self.getMatrix()).copy()
        tempMat[:, :, 0] = np.take(bCurve.mCurve, tempMat[:, :, 0])
        tempMat[:, :, 1] = np.take(gCurve.mCurve, tempMat[:, :, 1])
        tempMat[:, :, 2] = np.take(rCurve.mCurve, tempMat[:, :, 2])
        #Now we jimmy the np array into a cvMat
        image = cv.CreateImageHeader((tempMat.shape[1], tempMat.shape[0]), cv.IPL_DEPTH_8U, 3)
        cv.SetData(image, tempMat.tostring(), tempMat.dtype.itemsize * 3 * tempMat.shape[1])
        return Image(image, colorSpace=self._colorSpace)

    def applyIntensityCurve(self, curve):
        """
        Intensity applied to all three color channels

        Returns: Image
        """
        return self.applyRGBCurve(curve, curve, curve)

    def colorDistance(self, color = [0, 0, 0]):
        """
        Returns an image representing the distance of each pixel from a given color
        tuple, scaled between 0 (the given color) and 255.  Pixels distant from the
        given tuple will appear as brighter and pixels closest to the target color
        will be darker.

        By default this will give image intensity (distance from pure black)
        """
        pixels = np.array(self.getNumpy()).reshape(-1, 3)   #reshape our matrix to 1xN
        distances = spsd.cdist(pixels, [color]) #calculate the distance each pixel is
        distances *= (255.0/distances.max()) #normalize to 0 - 255
        return Image(distances.reshape(self.width, self.height)) #return an Image

    def erode(self, iterations=1):
        """
        Apply a morphological erosion. An erosion has the effect of removing small bits of noise
        and smothing blobs.
        This implementation uses the default openCV 3X3 square kernel
        Erosion is effectively a local minima detector, the kernel moves over the image and
        takes the minimum value inside the kernel.
        iterations - this parameters is the number of times to apply/reapply the operation
        See: http://en.wikipedia.org/wiki/Erosion_(morphology).
        See: http://opencv.willowgarage.com/documentation/cpp/image_filtering.html#cv-erode
        Example Use: A threshold/blob image has 'salt and pepper' noise.
        Example Code: ./examples/MorphologyExample.py

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        kern = cv.CreateStructuringElementEx(3, 3, 1, 1, cv.CV_SHAPE_RECT)
        cv.Erode(self.getBitmap(), retVal, kern, iterations)
        return Image(retVal, colorSpace=self._colorSpace)

    def dilate(self, iterations=1):
        """
        Apply a morphological dilation. An dilation has the effect of smoothing blobs while
        intensifying the amount of noise blobs.
        This implementation uses the default openCV 3X3 square kernel
        Erosion is effectively a local maxima detector, the kernel moves over the image and
        takes the maxima value inside the kernel.

        iterations - this parameters is the number of times to apply/reapply the operation

        See: http://en.wikipedia.org/wiki/Dilation_(morphology)
        See: http://opencv.willowgarage.com/documentation/cpp/image_filtering.html#cv-dilate
        Example Use: A part's blob needs to be smoother
        Example Code: ./examples/MorphologyExample.py

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        kern = cv.CreateStructuringElementEx(3, 3, 1, 1, cv.CV_SHAPE_RECT)
        cv.Dilate(self.getBitmap(), retVal, kern, iterations)
        return Image(retVal, colorSpace=self._colorSpace)

    def morphOpen(self):
        """
        morphologyOpen applies a morphological open operation which is effectively
        an erosion operation followed by a morphological dilation. This operation
        helps to 'break apart' or 'open' binary regions which are close together.

        See: http://en.wikipedia.org/wiki/Opening_(morphology)
        See: http://opencv.willowgarage.com/documentation/cpp/image_filtering.html#cv-morphologyex
        Example Use: two part blobs are 'sticking' together.
        Example Code: ./examples/MorphologyExample.py

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        temp = self.getEmpty()
        kern = cv.CreateStructuringElementEx(3, 3, 1, 1, cv.CV_SHAPE_RECT)
        cv.MorphologyEx(self.getBitmap(), retVal, temp, kern, cv.MORPH_OPEN, 1)
        return( Image(retVal) )

    def morphClose(self):
        """
        morphologyClose applies a morphological close operation which is effectively
        a dilation operation followed by a morphological erosion. This operation
        helps to 'bring together' or 'close' binary regions which are close together.

        See: http://en.wikipedia.org/wiki/Closing_(morphology)
        See: http://opencv.willowgarage.com/documentation/cpp/image_filtering.html#cv-morphologyex
        Example Use: Use when a part, which should be one blob is really two blobs.
        Example Code: ./examples/MorphologyExample.py

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        temp = self.getEmpty()
        kern = cv.CreateStructuringElementEx(3, 3, 1, 1, cv.CV_SHAPE_RECT)
        cv.MorphologyEx(self.getBitmap(), retVal, temp, kern, cv.MORPH_CLOSE, 1)
        return Image(retVal, colorSpace=self._colorSpace)

    def morphGradient(self):
        """
        The morphological gradient is the difference betwen the morphological
        dilation and the morphological gradient. This operation extracts the
        edges of a blobs in the image.

        See: http://en.wikipedia.org/wiki/Morphological_Gradient
        See: http://opencv.willowgarage.com/documentation/cpp/image_filtering.html#cv-morphologyex
        Example Use: Use when you have blobs but you really just want to know the blob edges.
        Example Code: ./examples/MorphologyExample.py

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        retVal = self.getEmpty()
        temp = self.getEmpty()
        kern = cv.CreateStructuringElementEx(3, 3, 1, 1, cv.CV_SHAPE_RECT)
        cv.MorphologyEx(self.getBitmap(), retVal, temp, kern, cv.MORPH_GRADIENT, 1)
        return Image(retVal, colorSpace=self._colorSpace )

    def histogram(self, numbins = 50):
        """
        Return a numpy array of the 1D histogram of intensity for pixels in the image
        Single parameter is how many "bins" to have.

        Returns: LIST
        """
        gray = self._getGrayscaleBitmap()

        (hist, bin_edges) = np.histogram(np.asarray(cv.GetMat(gray)), bins=numbins)
        return hist.tolist()

    def __getitem__(self, coord):
        ret = self.getMatrix()[tuple(reversed(coord))]
        if (type(ret) == cv.cvmat):
            (width, height) = cv.GetSize(ret)
            newmat = cv.CreateMat(height, width, ret.type)
            cv.Copy(ret, newmat) #this seems to be a bug in opencv
            #if you don't copy the matrix slice, when you convert to bmp you get
            #a slice-sized hunk starting at 0, 0
            return Image(newmat)
        return tuple(reversed(ret))

    def __setitem__(self, coord, value):
        value = tuple(reversed(value))  #RGB -> BGR
        if (is_tuple(self.getMatrix()[tuple(reversed(coord))])):
            self.getMatrix()[tuple(reversed(coord))] = value
        else:
            cv.Set(self.getMatrix()[tuple(reversed(coord))], value)
            self._clearBuffers("_matrix")

    def __sub__(self, other):
        newbitmap = self.getEmpty()
        if is_number(other):
            cv.SubS(self.getBitmap(), other, newbitmap)
        else:
            cv.Sub(self.getBitmap(), other.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __add__(self, other):
        newbitmap = self.getEmpty()
        if is_number(other):
            cv.AddS(self.getBitmap(), other, newbitmap)
        else:
            cv.Add(self.getBitmap(), other.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __and__(self, other):
        newbitmap = self.getEmpty()
        if is_number(other):
            cv.AndS(self.getBitmap(), other, newbitmap)
        else:
            cv.And(self.getBitmap(), other.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __or__(self, other):
        newbitmap = self.getEmpty()
        if is_number(other):
            cv.OrS(self.getBitmap(), other, newbitmap)
        else:
            cv.Or(self.getBitmap(), other.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __div__(self, other):
        newbitmap = self.getEmpty()
        if (not is_number(other)):
            cv.Div(self.getBitmap(), other.getBitmap(), newbitmap)
        else:
            cv.ConvertScale(self.getBitmap(), newbitmap, 1.0/float(other))
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __mul__(self, other):
        newbitmap = self.getEmpty()
        if (not is_number(other)):
            cv.Mul(self.getBitmap(), other.getBitmap(), newbitmap)
        else:
            cv.ConvertScale(self.getBitmap(), newbitmap, float(other))
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __pow__(self, other):
        newbitmap = self.getEmpty()
        cv.Pow(self.getBitmap(), newbitmap, other)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def __neg__(self):
        newbitmap = self.getEmpty()
        cv.Not(self.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def max(self, other):
        """
        The maximum value of my image, and the other image, in each channel
        If other is a number, returns the maximum of that and the number

        Returns: IMAGE
        """
        newbitmap = self.getEmpty()
        if is_number(other):
            cv.MaxS(self.getBitmap(), other.getBitmap(), newbitmap)
        else:
            cv.Max(self.getBitmap(), other.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def min(self, other):
        """
        The minimum value of my image, and the other image, in each channel
        If other is a number, returns the minimum of that and the number

        Returns: IMAGE
        """
        newbitmap = self.getEmpty()
        if is_number(other):
            cv.MaxS(self.getBitmap(), other.getBitmap(), newbitmap)
        else:
            cv.Max(self.getBitmap(), other.getBitmap(), newbitmap)
        return Image(newbitmap, colorSpace=self._colorSpace)

    def _clearBuffers(self, clearexcept = "_bitmap"):
        for k, v in self._initialized_buffers.items():
            if k == clearexcept:
                continue
            self.__dict__[k] = v

    def findBarcode(self, zxing_path = ""):
        """
        If you have the python-zxing library installed, you can find 2d and 1d
        barcodes in your image.  These are returned as Barcode feature objects
        in a FeatureSet.  The single parameter is the ZXing_path, if you
        don't have the ZXING_LIBRARY env parameter set.

        You can clone python-zxing at http://github.com/oostendo/python-zxing

        Returns: BARCODE
        """
        if not ZXING_ENABLED:
            return None

        if (not self._barcodeReader):
            if not zxing_path:
                self._barcodeReader = zxing.BarCodeReader()
            else:
                self._barcodeReader = zxing.BarCodeReader(zxing_path)

        tmp_filename = os.tmpnam() + ".png"
        self.save(tmp_filename)
        barcode = self._barcodeReader.decode(tmp_filename)
        os.unlink(tmp_filename)

        if barcode:
            return Barcode(self, barcode)
        else:
            return None

    #this function contains two functions -- the basic edge detection algorithm
    #and then a function to break the lines down given a threshold parameter
    def findLines(self, threshold=80, minlinelength=30, maxlinegap=10, cannyth1=50, cannyth2=100):
        """
        findLines will find line segments in your image and returns Line feature
        objects in a FeatureSet. The parameters are:
        * threshold, which determies the minimum "strength" of the line
        * min line length -- how many pixels long the line must be to be returned
        * max line gap -- how much gap is allowed between line segments to consider them the same line
        * cannyth1 and cannyth2 are thresholds used in the edge detection step, refer to _getEdgeMap() for details

        For more information, consult the cv.HoughLines2 documentation

        Returns: FEATURESET
        """
        em = self._getEdgeMap(cannyth1, cannyth2)

        lines = cv.HoughLines2(em, cv.CreateMemStorage(), cv.CV_HOUGH_PROBABILISTIC, 1.0, cv.CV_PI/180.0, threshold, minlinelength, maxlinegap)

        linesFS = FeatureSet()
        for l in lines:
            linesFS.append(Line(self, l))
        return linesFS

    def findChessboard(self, dimensions = (8, 5), subpixel = True):
        """
        Given an image, finds a chessboard within that image.  Returns the Chessboard featureset.
        The Chessboard is typically used for calibration because of its evenly spaced corners.

        The single parameter is the dimensions of the chessboard, typical one can be found in \SimpleCV\tools\CalibGrid.png

        returns a FeatureSet with the Chessboard feature, or none
        """
        corners = cv.FindChessboardCorners(self._getEqualizedGrayscaleBitmap(), dimensions, cv.CALIB_CB_ADAPTIVE_THRESH + cv.CALIB_CB_NORMALIZE_IMAGE + cv.CALIB_CB_FAST_CHECK )
        if(len(corners[1]) == dimensions[0]*dimensions[1]):
            if (subpixel):
                spCorners = cv.FindCornerSubPix(self.getGrayscaleMatrix(), corners[1], (11, 11), (-1, -1), (cv.CV_TERMCRIT_ITER | cv.CV_TERMCRIT_EPS, 10, 0.01))
            else:
                spCorners = corners[1]
            return FeatureSet([ Chessboard(self, dimensions, spCorners) ])
        else:
            return None

    def edges(self, t1=50, t2=100):
        """
        Finds an edge map Image using the Canny edge detection method.  Edges will be brighter than the surrounding area.

        The t1 parameter is roughly the "strength" of the edge required, and the value between t1 and t2 is used for edge linking.  For more information:

        <http://opencv.willowgarage.com/documentation/python/imgproc_feature_detection.html>
        <http://en.wikipedia.org/wiki/Canny_edge_detector>

        Returns: IMAGE
        """
        return Image(self._getEdgeMap(t1, t2), colorSpace=self._colorSpace)

    def _getEdgeMap(self, t1=50, t2=100):
        """
        Return the binary bitmap which shows where edges are in the image.  The two
        parameters determine how much change in the image determines an edge,
        and how edges are linked together.  For more information refer to:

        http://en.wikipedia.org/wiki/Canny_edge_detector
        http://opencv.willowgarage.com/documentation/python/imgproc_feature_detection.html?highlight=canny#Canny
        """

        if (self._edgeMap and self._cannyparam[0] == t1 and self._cannyparam[1] == t2):
            return self._edgeMap

        self._edgeMap = self.getEmpty(1)
        cv.Canny(self._getGrayscaleBitmap(), self._edgeMap, t1, t2)
        self._cannyparam = (t1, t2)

        return self._edgeMap

    def rotate(self, angle, mode="fixed", point=[-1, -1], scale = 1.0):
        """
        This rotates an image around a specific point by the given angle
        By default in "fixed" mode, the returned Image is the same dimensions as the original Image, and the contents will be scaled to fit.  In "full" mode the
        contents retain the original size, and the Image object will scale
        by default, the point is the center of the image.
        you can also specify a scaling parameter

        Returns: IMAGE
        """
        if( point[0] == -1 or point[1] == -1 ):
            point[0] = (self.width-1)/2
            point[1] = (self.height-1)/2

        if (mode == "fixed"):
            retVal = self.getEmpty()
            rotMat = cv.CreateMat(2, 3, cv.CV_32FC1)
            cv.GetRotationMatrix2D((float(point[0]), float(point[1])), float(angle), float(scale), rotMat)
            cv.WarpAffine(self.getBitmap(), retVal, rotMat)
            return Image(retVal, colorSpace=self._colorSpace)

        #otherwise, we're expanding the matrix to fit the image at original size
        rotMat = cv.CreateMat(2, 3, cv.CV_32FC1)
        # first we create what we thing the rotation matrix should be
        cv.GetRotationMatrix2D((float(point[0]), float(point[1])), float(angle), float(scale), rotMat)
        A = np.array([0, 0, 1])
        B = np.array([self.width, 0, 1])
        C = np.array([self.width, self.height, 1])
        D = np.array([0, self.height, 1])
        #So we have defined our image ABC in homogenous coordinates
        #and apply the rotation so we can figure out the image size
        a = np.dot(rotMat, A)
        b = np.dot(rotMat, B)
        c = np.dot(rotMat, C)
        d = np.dot(rotMat, D)
        #I am not sure about this but I think the a/b/c/d are transposed
        #now we calculate the extents of the rotated components.
        minY = min(a[1], b[1], c[1], d[1])
        minX = min(a[0], b[0], c[0], d[0])
        maxY = max(a[1], b[1], c[1], d[1])
        maxX = max(a[0], b[0], c[0], d[0])
        #from the extents we calculate the new size
        newWidth = np.ceil(maxX-minX)
        newHeight = np.ceil(maxY-minY)
        #now we calculate a new translation
        tX = 0
        tY = 0
        #calculate the translation that will get us centered in the new image
        if( minX < 0 ):
            tX = -1.0*minX
        elif(maxX > newWidth-1 ):
            tX = -1.0*(maxX-newWidth)

        if( minY < 0 ):
            tY = -1.0*minY
        elif(maxY > newHeight-1 ):
            tY = -1.0*(maxY-newHeight)

        #now we construct an affine map that will the rotation and scaling we want with the
        #the corners all lined up nicely with the output image.
        src = ((A[0], A[1]), (B[0], B[1]), (C[0], C[1]))
        dst = ((a[0]+tX, a[1]+tY), (b[0]+tX, b[1]+tY), (c[0]+tX, c[1]+tY))

        cv.GetAffineTransform(src, dst, rotMat)

        #calculate the translation of the corners to center the image
        #use these new corner positions as the input to cvGetAffineTransform
        retVal = cv.CreateImage((int(newWidth), int(newHeight)), 8, int(3))
        cv.WarpAffine(self.getBitmap(), retVal, rotMat)
        return Image(retVal, colorSpace=self._colorSpace)

    def rotate90(self):
        """
        Does a fast 90 degree rotation.
        """
        retVal = cv.CreateImage((self.height, self.width), cv.IPL_DEPTH_8U, 3)
        cv.Transpose(self.getBitmap(), retVal)
        return(Image(retVal, colorSpace=self._colorSpace))

    def shear(self, cornerpoints):
        """
        Given a set of new corner points in clockwise order, return a shear-ed Image
        that transforms the Image contents.  The returned image is the same
        dimensions.

        cornerpoints is a 2x4 array of point tuples

        Returns: IMAGE
        """
        src =  ((0, 0), (self.width-1, 0), (self.width-1, self.height-1))
        #set the original points
        aWarp = cv.CreateMat(2, 3, cv.CV_32FC1)
        #create the empty warp matrix
        cv.GetAffineTransform(src, cornerpoints, aWarp)

        return self.transformAffine(aWarp)

    def transformAffine(self, rotMatrix):
        """
        This helper function for shear performs an affine rotation using the supplied matrix.
        The matrix can be a either an openCV mat or an np.ndarray type.
        The matrix should be a 2x3

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        if(type(rotMatrix) == np.ndarray ):
            rotMatrix = npArray2cvMat(rotMatrix)
        cv.WarpAffine(self.getBitmap(), retVal, rotMatrix)
        return Image(retVal, colorSpace=self._colorSpace)

    def warp(self, cornerpoints):
        """
        Given a new set of corner points in clockwise order, return an Image with
        the images contents warped to the new coordinates.  The returned image
        will be the same size as the original image

        Returns: IMAGE
        """
        #original coordinates
        src = ((0, 0), (self.width-1, 0), (self.width-1, self.height-1), (0, self.height-1))

        pWarp = cv.CreateMat(3, 3, cv.CV_32FC1) #create an empty 3x3 matrix
        cv.GetPerspectiveTransform(src, cornerpoints, pWarp) #figure out the warp matrix

        return self.transformPerspective(pWarp)

    def transformPerspective(self, rotMatrix):

        """
        This helper function for warp performs an affine rotation using the supplied matrix.
        The matrix can be a either an openCV mat or an np.ndarray type.
        The matrix should be a 3x3

        Returns: IMAGE
        """
        retVal = self.getEmpty()
        if(type(rotMatrix) == np.ndarray ):
            rotMatrix = npArray2cvMat(rotMatrix)
        cv.WarpPerspective(self.getBitmap(), retVal, rotMatrix)
        return Image(retVal, colorSpace=self._colorSpace)

    def getPixel(self, x, y):
        """
        This function returns the RGB value for a particular image pixel given a specific row and column.
        """
        retVal = None
        if( x < 0 or x >= self.width ):
            warnings.warn("getRGBPixel: X value is not valid.")
        elif( y < 0 or y >= self.height ):
            warnings.warn("getRGBPixel: Y value is not valid.")
        else:
            retVal = cv.Get2D(self.getBitmap(), y, x)
        return retVal

    def getGrayPixel(self, x, y):
        """
        This function returns the Gray value for a particular image pixel given a specific row and column.
        """
        retVal = None
        if( x < 0 or x >= self.width ):
            warnings.warn("getGrayPixel: X value is not valid.")
        elif( y < 0 or y >= self.height ):
            warnings.warn("getGrayPixel: Y value is not valid.")
        else:
            retVal = cv.Get2D(self._getGrayscaleBitmap(), y, x)
            retVal = retVal[0]
        return retVal

    def getVertScanline(self, column):
        """
        This function returns a single column of RGB values from the image.
        """
        retVal = None
        if( column < 0 or column >= self.width ):
            warnings.warn("getVertRGBScanline: column value is not valid.")
        else:
            retVal = cv.GetCol(self.getBitmap(), column)
            retVal = np.array(retVal)
            retVal = retVal[:, 0, :]
        return retVal

    def getHorzScanline(self, row):
        """
        This function returns a single row of RGB values from the image.
        """
        retVal = None
        if( row < 0 or row >= self.height ):
            warnings.warn("getHorzRGBScanline: row value is not valid.")
        else:
            retVal = cv.GetRow(self.getBitmap(), row)
            retVal = np.array(retVal)
            retVal = retVal[0, :, :]
        return retVal

    def getVertScanlineGray(self, column):
        """
        This function returns a single column of gray values from the image.
        """
        retVal = None
        if( column < 0 or column >= self.width ):
            warnings.warn("getHorzRGBScanline: row value is not valid.")
        else:
            retVal = cv.GetCol(self._getGrayscaleBitmap(), column )
            retVal = np.array(retVal)
            #retVal = retVal.transpose()
        return retVal

    def getHorzScanlineGray(self, row):
        """
        This function returns a single row of RGB values from the image.
        """
        retVal = None
        if( row < 0 or row >= self.height ):
            warnings.warn("getHorzRGBScanline: row value is not valid.")
        else:
            retVal = cv.GetRow(self._getGrayscaleBitmap(), row )
            retVal = np.array(retVal)
            retVal = retVal.transpose()
        return retVal

    def crop(self, x , y, w, h, centered=False):
        """
        Crop attempts to use the x and y position variables and the w and h width
        and height variables to crop the image. When centered is false, x and y
        define the top and left of the cropped rectangle. When centered is true
        the function uses x and y as the centroid of the cropped region.

        The function returns a new image.
        """
        retVal = cv.CreateImage((w, h), cv.IPL_DEPTH_8U, 3)
        if( centered ):
            rectangle = (x-(w/2), y-(h/2), w, h)
        else:
            rectangle = (x, y, w, h)

        cv.SetImageROI(self.getBitmap(), rectangle)
        cv.Copy(self.getBitmap(), retVal)
        cv.ResetImageROI(self.getBitmap())
        return Image(retVal, colorSpace=self._colorSpace)

    def regionSelect(self, x1, y1, x2, y2 ):
        """
        Region select is similar to crop, but instead of taking a position and width
        and height values it simply takes to points on the image and returns the selected
        region. This is very helpful for creating interactive scripts that require
        the user to select a region.
        """
        w = abs(x1-x2)
        h = abs(y1-y2)

        retVal = None
        if( w <= 0 or h <= 0 or w > self.width or h > self.height ):
            warnings.warn("regionSelect: the given values will not fit in the image or are too small.")
        else:
            xf = x2
            if( x1 < x2 ):
                xf = x1
            yf = y2
            if( y1 < y2 ):
                yf = y1
            retVal = self.crop(xf, yf, w, h)

        return retVal

    def clear(self):
        """
        This is a slightly unsafe method that clears out the entire image state
        it is usually used in conjunction with the drawing blobs to fill in draw
        only a single large blob in the image.
        """
        cv.SetZero(self._bitmap)
        self._clearBuffers()

    def drawText(self, text = "", x = None, y = None, color = Color.BLUE, fontsize = 16):
        """
        This function draws the string that is passed on the screen at the specified coordinates

        The Default Color is blue but you can pass it various colors
        The text will default to the center of the screen if you don't pass it a value

        returns Image

        """
        if(x == None):
            x = (self.width / 2)
        if(y == None):
            y = (self.height / 2)

        self.getDrawingLayer().setFontSize(fontsize)
        self.getDrawingLayer().text(text, (x, y), color)

    def show(self, type = 'window'):
        """
        This function automatically pops up a window and shows the current image
        """
        if(type == 'browser'):
          import webbrowser
          js = JpegStreamer(8080)
          self.save(js)
          webbrowser.open("http://localhost:8080", 2)
          return js
        elif (type == 'window'):
          from SimpleCV.Display import Display
          d = Display(self.size())
          self.save(d)
          return d
        else:
          print "Unknown type to show"

    def _surface2Image(self,surface):
        imgarray = pg.surfarray.array3d(surface)
        retVal = Image(imgarray)
        retVal._colorSpace = ColorSpace.RGB
        return retVal.toBGR().rotate90()

    def _image2Surface(self,img):
        return pg.image.fromstring(img.getPIL().tostring(),img.size(), "RGB")
        #return pg.surfarray.make_surface(img.toRGB().getNumpy())

    def toPygameSurface(self):
        """
        Converts this image to a pygame surface. This is useful if you want
        to treat an image as a sprite to render onto an image. An example
        would be rendering blobs on to an image. THIS IS EXPERIMENTAL.
        """
        return pg.image.fromstring(self.getPIL().tostring(),self.size(), "RGB")

    def addDrawingLayer(self, layer = ""):
        """
        Push a new drawing layer onto the back of the layer stack
        """
        if not layer:
            layer = DrawingLayer(self.size())
        self._mLayers.append(layer)
        return len(self._mLayers)-1

    def insertDrawingLayer(self, layer, index):
        """
        Insert a new layer into the layer stack at the specified index
        """
        self._mLayers.insert(index, layer)
        return None

    def removeDrawingLayer(self, index):
        """
        Remove a layer from the layer stack based on the layer's index.
        """
        return self._mLayers.pop(index)

    def getDrawingLayer(self, index = -1):
        """
        Return a drawing layer based on the provided index.  If not provided, will
        default to the top layer.  If no layers exist, one will be created
        """
        if not len(self._mLayers):
            self.addDrawingLayer()

        return self._mLayers[index]

    def dl(self, index = -1):
        """
        Alias for getDrawingLayer()
        """
        return self.getDrawingLayer(index)

    def clearLayers(self):
        """
        Remove all of the drawing layers.
        """
        for i in self._mLayers:
            self._mLayers.remove(i)

        return None

        #render the image.
    def _renderImage(self, layer):
        imgSurf = self.getPGSurface(self).copy()
        imgSurf.blit(layer._mSurface, (0, 0))
        return Image(imgSurf)

    def applyLayers(self, indicies=-1):
        """
        Render all of the layers onto the current image and return the result.
        Indicies can be a list of integers specifying the layers to be used.
        """
        if not len(self._mLayers):
            return self

        final = DrawingLayer((self.width, self.height))
        if(indicies==-1 and len(self._mLayers) > 0 ):
            #retVal = self
            self._mLayers.reverse()
            for layers in self._mLayers: #compose all the layers
                layers.renderToOtherLayer(final)
            self._mLayers.reverse()
            #then draw them
            imgSurf = self.getPGSurface().copy()
            imgSurf.blit(final._mSurface, (0, 0))
            return Image(imgSurf)
        else:
            retVal = self
            indicies.reverse()
            for idx in indicies:
                retVal = self._mLayers[idx].renderToOtherLayer(final)
            imgSurf = self.getPGSurface().copy()
            imgSurf.blit(final._mSurface, (0, 0))
            indicies.reverse()
            return Image(imgSurf)

from SimpleCV.BlobMaker import BlobMaker

max_stars_repo_name: quantombone/SimpleCV

max_issues_repo_name: quantombone/SimpleCV

max_forks_repo_name: quantombone/SimpleCV


text:
# %% [markdown]
# # The Cancer Gene Atlas (TCGA)
# %% [markdown]
# Created by: Charles Dai <br>
# Credit to: Moshe Silverstein
#
# Data Source: https://gdc.cancer.gov/
# %%
# appyter init
from appyter import magic
magic.init(lambda _=globals: _())
# %%
import sys
import os
from datetime import date
import gzip
import io
import requests
import json

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
from tqdm import tqdm

import harmonizome.utility_functions as uf
import harmonizome.lookup as lookup
# %%
# from clustergrammer_widget import *
# net = Network(clustergrammer_widget)
# %%
%load_ext autoreload
%autoreload 2
# %% [markdown]
# ### Notebook Information
# %%
print('This notebook was run on:', date.today(), '\nPython version:', sys.version)
# %% [markdown]
# # Initialization
# %%
%%appyter hide_code
{% do SectionField(
    name='data',
    title='Load Data',
    subtitle='Upload Files from the Comparative Toxicogenomics Database',
) %}
project_name = 'Mixed germ cell tumor'
# %%
base_url = 'https://api.gdc.cancer.gov/'
files_endpt = base_url + 'files/'
cases_endpt = base_url + 'cases/'
# %% [markdown]
# ### Load Mapping Dictionaries
# %%
symbol_lookup, geneid_lookup = lookup.get_lookups()
# %% [markdown]
# ### Output Path
# %%
output_name = 'tcga'

path = 'Output/TCGA'
if not os.path.exists(path):
    os.makedirs(path)
# %% [markdown]
# # Load Data
# %% [markdown]
# ## File Name and IDs
# %%
fields = [
    'cases.case_id'
]

filters = {
    'op': 'and',
    'content': [{
        'op': 'in',
        'content': {
            'field': 'experimental_strategy',
            'value': ['RNA-Seq'],
        }
    },
    {
        'op': 'in',
        'content': {
            'field': 'access',
            'value': ['open'],
        }
    },
    {
        'op': 'in',
        'content': {
            'field': 'file_name',
            'value': ['*htseq.counts.gz'],
        }
    },
    {
        'op': 'in',
        'content': {
            'field': 'cases.project.name',
            'value': [project_name],
        }
    }
    ],
}
# %%
params = {
    'fields': ','.join(fields),
    'filters': json.dumps(filters),
    'size': 100000,
    'facets': 'cases.case_id'
}
response = requests.get('https://api.gdc.cancer.gov/files', params=params)
data = response.json()['data']['hits']
# %%
files = pd.DataFrame([(f['id'], f['cases'][0]['case_id']) for f in data], columns=['file_id', 'case_id']).set_index('file_id')
print(files.shape)
files.head()
# %%
matrix = pd.DataFrame()

for file_id in tqdm(files.index, unit='samples'):
    response = requests.get(data_endpt + file_id, headers = {"Content-Type": "application/json"})
    string_data = io.StringIO(str(gzip.decompress(response.content), 'utf-8'))
    matrix = pd.concat([matrix, pd.read_csv(string_data, sep='\t', header=None, names=['ENSMBL ID', files.loc[file_id, 'case_id']], index_col=0)], axis=1)
matrix.head()
# %%
matrix.index = matrix.index.map(lambda x: x.split('.')[0])
matrix.head()
# %%
cases_fields = requests.get(cases_endpt + '_mapping').json()['fields']
keyfields = [field for field in cases_fields if
    any(word in field for word in ['demographic', 'diagnoses']) and 'treatment' not in field]
# %%
sample_meta = pd.DataFrame()
for case_id in tqdm(files['case_id'].drop_duplicates(), unit='cases'):
    response = requests.get(cases_endpt + case_id, params={'fields': ','.join(keyfields)})
    data = response.json()['data']
    sample = pd.DataFrame([{'case_id': case_id, **data['demographic'], **data['diagnoses'][0]}])
    sample_meta = pd.concat([sample_meta, sample])
sample_meta = sample_meta.set_index('case_id')
sample_meta.head()
# %% [markdown]
# # Pre-process Data
# %% [markdown]
# ## Create Matrix of Inference Scores
# %%
# Remove duplicates by averaging first
matrix = matrix.set_index(['# GeneSymbol', 'DiseaseName'])
matrix = matrix.groupby(level=[0,1]).mean().unstack()
matrix.head()
# %%
matrix.columns = matrix.columns.droplevel(0)
matrix.columns.name = 'Disease Name'
matrix.index.name = 'Gene Symbol'
matrix.head()
# %%
matrix.shape
# %% [markdown]
# ## Save Unfiltered Matrix to file
# %%
uf.save_data(matrix, path, output_name + '_matrix_unfiltered',
            compression='gzip', dtype=np.float32)
# %% [markdown]
# # Filter Data
# %% [markdown]
# ## Map Gene Symbols to Up-to-date Approved Gene Symbols
# %%
matrix = uf.map_symbols(matrix, symbol_lookup)
matrix.shape
# %% [markdown]
# ## Merge Duplicate Genes By Rows and Duplicate Columns
# %%
matrix = uf.merge(matrix, 'row')
matrix = uf.merge(matrix, 'column')
matrix.shape
# %% [markdown]
# ## Remove Genes that are More Than 95% Missing or Zero Inference
# %%
matrix = matrix.replace(0, np.nan).dropna(
    thresh=0.05 * matrix.shape[1], axis=0).replace(np.nan, 0)
matrix.head()
# %%
matrix.shape
# %% [markdown]
# ## Normalize Matrix (Quantile Normalize the Matrix by Column)
# %%
matrix = uf.quantile_normalize(matrix)
matrix.head()
# %% [markdown]
# ## Normalize Matrix (Z-Score the Rows)
# %%
matrix = uf.zscore(matrix)
matrix.head()
# %% [markdown]
# ## Histogram of First Sample
# %%
matrix.iloc[:, 0].hist(bins=100)
# %% [markdown]
# ## Histogram of First Gene
# %%
matrix.iloc[0, :].hist(bins=100)
# %% [markdown]
# ## Save Filtered Matrix
# %%
uf.save_data(matrix, path, output_name + '_matrix_filtered',
            ext='tsv', compression='gzip')
# %% [markdown]
# # Analyze Data
# %% [markdown]
# ## Create Gene List
# %%
gene_list = uf.gene_list(matrix, geneid_lookup)
gene_list.head()
# %%
gene_list.shape
# %%
uf.save_data(gene_list, path, output_name + '_gene_list',
            ext='tsv', compression='gzip', index=False)
# %% [markdown]
# ## Create Attribute List
# %%
attribute_list = uf.attribute_list(matrix)
attribute_list.head()
# %%
attribute_list.shape
# %%
uf.save_data(attribute_list, path, output_name + '_attribute_list',
            ext='tsv', compression='gzip')
# %% [markdown]
# ## Create matrix of Standardized values (values between -1, and 1)
# %%
standard_matrix = uf.standardized_matrix(matrix)
standard_matrix.head()
# %%
uf.save_data(standard_matrix, path, output_name + '_standard_matrix',
            ext='tsv', compression='gzip')
# %% [markdown]
# ## Plot of A Single Celltype, Normalized Value vs. Standardized Value
# %%
plt.plot(matrix[matrix.columns[0]],
         standard_matrix[standard_matrix.columns[0]], 'bo')
plt.xlabel('Normalized Values')
plt.ylabel('Standardized Values')
plt.title(standard_matrix.columns[0])
plt.grid(True)
# %% [markdown]
# ## Create Ternary Matrix
# %%
ternary_matrix = uf.ternary_matrix(standard_matrix)
ternary_matrix.head()
# %%
uf.save_data(ternary_matrix, path, output_name + '_ternary_matrix',
            ext='tsv', compression='gzip')
# %% [markdown]
# ## Create Gene and Attribute Set Libraries
# %%
uf.save_setlib(ternary_matrix, 'gene', 'up', path, output_name + '_gene_up_set')
# %%
uf.save_setlib(ternary_matrix, 'gene', 'down', path, output_name + '_gene_down_set')
# %%
uf.save_setlib(ternary_matrix, 'attribute', 'up', path,
                           output_name + '_attribute_up_set')
# %%
uf.save_setlib(ternary_matrix, 'attribute', 'down', path,
                             output_name + '_attribute_down_set')
# %% [markdown]
# ## Create Attribute Similarity Matrix
# %%
attribute_similarity_matrix = uf.similarity_matrix(standard_matrix.T, 'cosine')
attribute_similarity_matrix.head()
# %%
uf.save_data(attribute_similarity_matrix, path,
            output_name + '_attribute_similarity_matrix',
            compression='npz', symmetric=True, dtype=np.float32)
# %%
# net.load_df(attribute_similarity_matrix.iloc[:,:].copy())
# net.filter_N_top('row', rank_type='sum', N_top=300)
# net.cluster()
# net.widget()
# %% [markdown]
# ## Create Gene Similarity Matrix
# %%
gene_similarity_matrix = uf.similarity_matrix(standard_matrix, 'cosine')
gene_similarity_matrix.head()
# %%
uf.save_data(gene_similarity_matrix, path,
            output_name + '_gene_similarity_matrix',
            compression='npz', symmetric=True, dtype=np.float32)
# %% [markdown]
# ## Create Gene-Attribute Edge List
# %%
edge_list = uf.edge_list(standard_matrix)
uf.save_data(edge_list, path, output_name + '_edge_list',
        ext='tsv', compression='gzip')
# %% [markdown]
# # Create Downloadable Save File
# %%
uf.archive(path)
# %% [markdown]
# ### Link to download output files: [click here](./output_archive.zip)

ext: py

reason: import numpy


text:
# encoding=utf-8

import os
import numpy as np
from Moldata import Moldata
from Reaction import Reaction
import matplotlib.pyplot as plt

def shermo(
    input,
    Temp=300.0,
    sclZPE=1.0,
    sclHeat=1.0,
    sclS=1.0,
    E='N/A',
    shermo_path='Shermo'
    ):

    if E=='N/A':
        command=shermo_path+' '+input+' -T '+str(Temp)+' -P 1.0 -sclZPE '+str(sclZPE)+' -sclheat '+str(sclHeat)+' -sclS '+str(sclS)
    else:
        command=shermo_path+' '+input+' -E '+str(E)+' -T '+str(Temp)+' -P 1.0 -sclZPE '+str(sclZPE)+' -sclheat '+str(sclHeat)+' -sclS '+str(sclS)

    os.system(command+'> tmp')

    file=open('tmp','r')
    tmp=file.readlines()[-13:]
    file.close()
    os.remove('tmp')

    data=['','','']
    data[0]=tmp[-4].split()
    data[1]=tmp[-1].split()
    data[2]=tmp[0].split()

    result={}
    result['U0K']=data[0][-2]
    result['GTK']=data[1][-2]
    result['Q']=data[2][-1]
    return result

def run_shermo(
    Nmol,
    Rinput='',
    R2input='N/A',
    TSinput='',
    Pinput='',
    P2input='N/A',
    iFreq=-0.0,
    Temp=300.0,
    RE='N/A',
    R2E='N/A',
    TSE='N/A',
    PE='N/A',
    P2E='N/A',
    sclZPE=1.0,
    sclHeat=1.0,
    sclS=1.0,
    shermo_path='Shermo'
    ):

    if Nmol==1:
        Rresult=shermo(
            input=Rinput,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS, E=RE,
            shermo_path=shermo_path
            )

        TSresult=shermo(
            input=TSinput,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            E=TSE,
            shermo_path=shermo_path
            )

        Presult=shermo(
            input=Pinput,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            E=PE,
            shermo_path=shermo_path
            )

        R=Moldata(
            U0K=Rresult['U0K'],
            GTK=Rresult['GTK'],
            Q=Rresult['Q'],
            EUnit='Eh'
            )

        TS=Moldata(
            U0K=TSresult['U0K'],
            GTK=TSresult['GTK'],
            Q=TSresult['Q'],
            EUnit='Eh'
            )

        if P2input=='N/A':
            P=Moldata(
                U0K=Presult['U0K'],
                GTK=Presult['GTK'],
                Q=Presult['Q'],
                EUnit='Eh'
                )
        else:
            P2result=shermo(
                input=Pinput,
                Temp=Temp,
                sclZPE=sclZPE,
                sclHeat=sclHeat,
                sclS=sclS,
                E=PE,
                shermo_path=shermo_path
                )

            P1=Moldata(
                U0K=Presult['U0K'],
                GTK=Presult['GTK'],
                Q=Presult['Q'],
                EUnit='Eh'
                )

            P2=Moldata(
                U0K=P2result['U0K'],
                GTK=P2result['GTK'],
                Q=P2result['Q'],
                EUnit='Eh'
                )

            P=P1+P2

        cal=Reaction(
            Nmol=1,
            molR=R,
            molTS=TS,
            molP=P,
            Temp=Temp,
            iFreq=iFreq
            )

        cal.printf()

    elif Nmol==2 and R2input!='N/A':
        Rresult=shermo(
            input=Rinput,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            E=RE,
            shermo_path=shermo_path
            )

        R2result=shermo(
            input=R2input,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            E=R2E,
            shermo_path=shermo_path
            )

        TSresult=shermo(
            input=TSinput,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            E=TSE,
            shermo_path=shermo_path
            )

        Presult=shermo(
            input=Pinput,
            Temp=Temp,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            E=PE,
            shermo_path=shermo_path
            )

        R=Moldata(
            U0K=Rresult['U0K'],
            GTK=Rresult['GTK'],
            Q=Rresult['Q'],
            EUnit='Eh'
            )

        R2=Moldata(
            U0K=R2result['U0K'],
            GTK=R2result['GTK'],
            Q=R2result['Q'],
            EUnit='Eh'
            )

        TS=Moldata(
            U0K=TSresult['U0K'],
            GTK=TSresult['GTK'],
            Q=TSresult['Q'],
            EUnit='Eh'
            )

        if P2input=='N/A':
            P=Moldata(
                U0K=Presult['U0K'],
                GTK=Presult['GTK'],
                Q=Presult['Q'],
                EUnit='Eh'
                )
        else:
            P2result=shermo(
                input=Pinput,
                Temp=Temp,
                sclZPE=sclZPE,
                sclHeat=sclHeat,
                sclS=sclS,
                E=PE,
                shermo_path=shermo_path
                )

            P1=Moldata(
                U0K=Presult['U0K'],
                GTK=Presult['GTK'],
                Q=Presult['Q'],
                EUnit='Eh'
                )

            P2=Moldata(
                U0K=P2result['U0K'],
                GTK=P2result['GTK'],
                Q=P2result['Q'],
                EUnit='Eh'
                )

            P=P1+P2

        cal=Reaction(
            Nmol=2,
            molR=R+R2,
            molTS=TS,
            molP=P,
            Temp=Temp,
            iFreq=iFreq
            )

        cal.printf()

    elif Nmol==2 and R2input=='N/A':
        print("Parameter R2input is necessary! ")
        exit()
    else:
        print("The number of molecules in primitive chemical reaction is not suitable! ")
        exit()
    return cal

def scan_shermo(
    Nmol,
    Rinput='',
    R2input='N/A',
    TSinput='',
    Pinput='',
    P2input='N/A',
    iFreq=-0.0,
    Temp=[],
    RE='N/A',
    R2E='N/A',
    TSE='N/A',
    PE='N/A',
    P2E='N/A',
    sclZPE=1.0,
    sclHeat=1.0,
    sclS=1.0,
    shermo_path='Shermo'
    ):

    result=[]
    Xaxis=Temp
    GYaxis=[]
    QYaxis=[]

    for temp in Temp:
        calc=run_shermo(
            Nmol=Nmol,
            Rinput=Rinput,
            R2input=R2input,
            TSinput=TSinput,
            Pinput=Pinput,
            P2input=P2input,
            iFreq=iFreq,
            Temp=float(temp),
            RE=RE,
            R2E=R2E,
            TSE=TSE,
            PE=PE,
            P2E=P2E,
            sclZPE=sclZPE,
            sclHeat=sclHeat,
            sclS=sclS,
            shermo_path=shermo_path
            )
        GYaxis.append(calc.get_kinetic_g())
        QYaxis.append(calc.get_kinetic_q())

    xx=np.linspace(min(Temp),max(Temp),200)
    GFunc=np.poly1d(np.polyfit(Xaxis,GYaxis,3))
    QFunc=np.poly1d(np.polyfit(Xaxis,QYaxis,3))

    result.append(Temp)
    result.append(GYaxis)
    result.append(QYaxis)

    if Nmol==1:
        print(f'''

---------------------------------------------------------
                        Part IV
  Relationship between Reaction Kinetics and Temperature
---------------------------------------------------------

Gibbs Free Energy Method:

    k(T)/s-1

    =  ({GFunc.coeffs[0]:.3E})\t * (T/K)^3

     + ({GFunc.coeffs[1]:.3E})\t * (T/K)^2

     + ({GFunc.coeffs[2]:.3E})\t * (T/K)

     + ({GFunc.coeffs[3]:.3E})

    RMSD=

Partition Function Method:

    k(T)/s-1

    =  ({QFunc.coeffs[0]:.3E})\t * (T/K)^3

     + ({QFunc.coeffs[1]:.3E})\t * (T/K)^2

     + ({QFunc.coeffs[2]:.3E})\t * (T/K)

     + ({QFunc.coeffs[3]:.3E})

    RMSD=

---------------------------------------------------------
        ''')

    elif Nmol==2:
        print(f'''

---------------------------------------------------------
                        Part IV
  Relationship between Reaction Kinetics and Temperature
---------------------------------------------------------

Gibbs Free Energy Method:

    k(T)/((mol/L)-1*s-1)

    =  ({GFunc.coeffs[0]:.3E})\t * (T/K)^3

     + ({GFunc.coeffs[1]:.3E})\t * (T/K)^2

     + ({GFunc.coeffs[2]:.3E})\t * (T/K)

     + ({GFunc.coeffs[3]:.3E})

    RMSD=

Partition Function Method:

    k(T)/((mol/L)-1*s-1)

    =  ({QFunc.coeffs[0]:.3E})\t * (T/K)^3

     + ({QFunc.coeffs[1]:.3E})\t * (T/K)^2

     + ({QFunc.coeffs[2]:.3E})\t * (T/K)

     + ({QFunc.coeffs[3]:.3E})

    RMSD=

---------------------------------------------------------
        ''')

    else:
        print("The number of molecules in primitive chemical reaction is not suitable! ")
        exit()

    print('''

d8b   db  .d88b.  d8888b. .88b  d88.  .d8b.  db
888o  88 .8P  Y8. 88  `8D 88'YbdP`88 d8' `8b 88
88V8o 88 88    88 88oobY' 88  88  88 88ooo88 88
88 V8o88 88    88 88`8b   88  88  88 88~~~88 88
88  V888 `8b  d8' 88 `88. 88  88  88 88   88 88booo.
VP   V8P  `Y88P'  88   YD YP  YP  YP YP   YP Y88888P

d88888b d8b   db d8888b.      db
88'     888o  88 88  `8D      88
88ooooo 88V8o 88 88   88      YP
88~~~~~ 88 V8o88 88   88
88.     88  V888 88  .8D      db
Y88888P VP   V8P Y8888D'      YP

            ''')

    plt.plot(Xaxis,GYaxis,'rx',xx,GFunc(xx),'-')
    plt.title('Reaction Rate Constant (Gibbs Free Energy Method)')
    if Nmol==1:
        plt.ylabel('k/(s-1)')
    elif Nmol==2:
        plt.ylabel('k/((mol/L)-1*s-1)')
    else:
        print("The number of molecules in primitive chemical reaction is not suitable! ")
        exit()
    plt.xlabel('Temperature/K')
    plt.figure(2)
    plt.plot(Xaxis,QYaxis,'rx',xx,QFunc(xx),'-')
    plt.title('Reaction Rate Constant (Partition Function Method)')
    if Nmol==1:
        plt.ylabel('k/(s-1)')
    elif Nmol==2:
        plt.ylabel('k/((mol/L)-1*s-1)')
    else:
        print("The number of molecules in primitive chemical reaction is not suitable! ")
        exit()
    plt.xlabel('Temperature/K')
    plt.show()

    return result

if __name__ == "__main__":
    result=scan_shermo(
        Nmol=2,
        Rinput='DA-R.log',
        R2input='DA-R2.log',
        TSinput='DA-TS.log',
        Pinput='DA-P.log',
        iFreq=-2000.2842,
        Temp=list(range(273,373))
        )
    print(result)

max_stars_repo_name: Linqiaosong/PyQTST

max_issues_repo_name: Linqiaosong/PyTST

max_forks_repo_name: Linqiaosong/PyTST


text:
import random
import json
import numpy as np
import spacy
import ast

from .models import Scotch, EmotionSyn, EmotionScotchSimilarity, Summaries

class RecommendationEngine():

	def match_emotion(self, emotion):
		# dict to store similarity values
		similarity = {}

		for e in EmotionSyn.objects.all():
			# check and see if emotion matches one in dictionary
			if emotion == e.parent_emotion or emotion in ast.literal_eval(e.synonyms):
				# if so, return emotion
				return e.parent_emotion
			else:
				pass

		# if not, calculate similarity between emotion and dictionary
		for e in EmotionSyn.objects.all():
			# initialize English model
			nlp = spacy.load('en_core_web_sm')

			doc1 = nlp(emotion)
			doc2 = nlp(e.parent_emotion)
			similarity[e.parent_emotion] = doc1.similarity(doc2)

			# return max
			return max(similarity, key=similarity.get)

	def get_rec(self, emotion):

		# get the emotion for recommendation
		self.rec_emotion = self.match_emotion(emotion)
		# get index of emotion
		self.rec_idx = EmotionSyn.objects.get(parent_emotion=self.rec_emotion).pk
		# extract similarities as list
		similarities = ast.literal_eval(EmotionScotchSimilarity.objects.get(pk=self.rec_idx).similarity)
		# convert to numpy array
		sim_array = np.asarray(similarities)
		# get indices of top three scotches based on precomputed similarity scores
		scotch_idx = random.sample(list(sim_array.argsort()[-25:][::-1]), 1)[0]

		return scotch_idx


text:
import os
import sys

# Fisrt, we add the location of the library to test to the PYTHON path
if sys.platform=="cygwin":
    from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
    from liblibra_core import *

from libra_py import hpc_utils
from libra_py import data_read
from libra_py import data_outs
from libra_py import units
from libra_py import QE_methods
from libra_py.workflows.nbra import step3
#import libra_py.workflows.nbra.step2_analysis as step2a

#import numpy as np
#from matplotlib.mlab import griddata

import matplotlib.pyplot as plt   # plots
#%matplotlib inline

plt.rc('axes', titlesize=24)      # fontsize of the axes title
plt.rc('axes', labelsize=20)      # fontsize of the x and y labels
plt.rc('legend', fontsize=20)     # legend fontsize
plt.rc('xtick', labelsize=16)    # fontsize of the tick labels
plt.rc('ytick', labelsize=16)    # fontsize of the tick labels

plt.rc('figure.subplot', left=0.2)
plt.rc('figure.subplot', right=0.95)
plt.rc('figure.subplot', bottom=0.13)
plt.rc('figure.subplot', top=0.88)

colors = {}

colors.update({"11": "#8b1a0e"})  # red
colors.update({"12": "#FF4500"})  # orangered
colors.update({"13": "#B22222"})  # firebrick
colors.update({"14": "#DC143C"})  # crimson

colors.update({"21": "#5e9c36"})  # green
colors.update({"22": "#006400"})  # darkgreen
colors.update({"23": "#228B22"})  # forestgreen
colors.update({"24": "#808000"})  # olive

colors.update({"31": "#8A2BE2"})  # blueviolet
colors.update({"32": "#00008B"})  # darkblue

colors.update({"41": "#2F4F4F"})  # darkslategray

clrs_index = ["11", "21", "31", "41", "12", "22", "32", "13","23", "14", "24"]

"""
  Total number of electrons 58 => 28 alp + 28 bet

                                            (HOMO)  (LUMO)
  1, 2, .. [ 20, 21, 22, 23, 24, 25, 26, 27, 28,    29   , ....    39,] 40, 41, ....

  Output matrices:

  [20, 21, ... 27, 28, 29, 30, ...  39, ] [40, 41, ... 48,         ... 59 ]

   0    1     [ 7,  8,  9,  10,]      19   20  21   [27,28, 29, 30,]  ... 39

                1   2   3    4                      -5. -6. -7. -8

"""

params = { "data_set_paths" : [os.getcwd()+"/traj1_in/"],
           "data_dim":40, "active_space":range(0,40),
           "isnap":0,  "fsnap":29,
           "data_re_prefix" : "S_dia_ks_", "data_re_suffix" : "_re",
           "data_im_prefix" : "S_dia_ks_", "data_im_suffix" : "_im"  }
S = data_read.get_data_sets(params)

params.update({ "data_re_prefix" : "St_dia_ks_", "data_re_suffix" : "_re",
                "data_im_prefix" : "St_dia_ks_", "data_im_suffix" : "_im"  } )
St = data_read.get_data_sets(params)

params.update({ "data_re_prefix" : "hvib_dia_", "data_re_suffix" : "_re",
                "data_im_prefix" : "hvib_dia_", "data_im_suffix" : "_im"  } )
Hvib_ks = data_read.get_data_sets(params)

# Remove the previous results and temporary working directory from the previous runs
os.system("rm -r traj1_out")

# Create the new results directory
os.system("mkdir traj1_out")

params = { "SD_basis" : [ [ 7,  8, -27, -28], [7, 8, -27, -29], [7, 8, -27, -30] ],
           "SD_energy_corr" : [0.0, 0.0, 0.0],
           "CI_basis" : [ [1.0, 0.0, 0.0],
                          [0.0, 1.0, 0.0],
                          [0.0, 0.0, 1.0]
                        ],
          "output_set_paths" : [os.getcwd()+"/traj1_out/"],
          "dt" : 1.0*units.fs2au,
          "do_orthogonalization":1,
          "do_state_reordering":2,
          "state_reordering_alpha":0.0,
          "do_phase_correction":1,
          "do_output" : 1,
          "Hvib_re_prefix":"Hvib_",  "Hvib_re_suffix":"_re",
          "Hvib_im_prefix":"Hvib_",  "Hvib_im_suffix":"_im"
         }

Hvib = step3.run(S, St, Hvib_ks, params)


text:
import numpy as np
import pandas as pd

class make_fisher_info_matrix():
  def __init__(self, dataDir, testCosts):
    self.txset = pd.read_csv(f"{ dataDir }/ratrtpcrabvec.csv")
    self.highs = pd.read_csv(f"{ dataDir }/ratrtpcrabhighs.csv")
    self.highs = self.highs.drop(columns=['Infection.State']).to_numpy()
    self.testCosts = testCosts

    self.Spec = [0.975,0.97,0.977] # Spec is the specificity values for  0=RAT,1=RTPCR,2=AB
    self.Sen = [0.5,0.95,0.921]   # Sen are the sensitivity values for  0=RAT,1=RTPCR,2=AB

    self.ttypevalid, self.patternCost = self.getCosting()
    self.patternCost = self.patternCost[1:]

    self.lenw = len(self.patternCost)

  def set_sensitivities(self, senRAT, senRTPCR, senIGG):
    self.Sen = [senRAT, senRTPCR, senIGG]

  def set_specificities(self, specRAT, specRTPCR, specIGG):
    self.Spec = [specRAT, specRTPCR, specIGG]

  def get_lenw(self):
    return self.lenw

  def get_patternCost(self):
    return self.patternCost

  def get_ttypevalid(self):
    return self.ttypevalid

  def getCosting(self):
    ttypevalid = self.txset.copy()

    #null value based values
    ttypevalid['RAT'] = 1 - ttypevalid['RAT'].apply(lambda x: 1 if pd.isna(x) else 0)
    ttypevalid['RTPCR'] = 1 - ttypevalid['RTPCR'].apply(lambda x: 1 if pd.isna(x) else 0)
    ttypevalid['AB'] = 1 - ttypevalid['AB'].apply(lambda x: 1 if pd.isna(x) else 0)

    ttypevalid = ttypevalid.drop_duplicates() # get unique rows
    ttypevalid = ttypevalid.drop(columns=['Test.Type']).to_numpy() #make a 2D matrix

    return ttypevalid, np.matmul(ttypevalid, np.array(self.testCosts)) #matrix multiplilcation

  # # This function estimates q_t(x;p), the probability of seeing test outcome x
  # # when test type is t and the disease prevalence vector is p.
  # # p = c(p[1], p[2], p[3]). ActiveInfection, ActiveABonly, ActiveForBoth
  # # Last state carries the remaining probability.
  def qxp(self, x, p):
      # Structure of x and p : x = ("0", "1", ""), p = (0.1, 0.2, 0.3)x
      prod = [1, 1, 1, 1]

      for j in range(3):
          if(x[j] == "0.0"):
              for i in range(4):
                  # For state i, test j, if correct response highs[i,j] is 1,
                  # this is a sensitivity error.
                  # If correct response is 0, this is correct specificity.
                  prod[i] = (prod[i])*(self.highs[i, j]*(1-self.Sen[j]) + (1-self.highs[i, j])*self.Spec[j])

          if(x[j] == "1.0"):
              for i in range(4):
                  # See explanation above
                  prod[i] = prod[i]*(self.highs[i, j]*(self.Sen[j]) + (1-self.highs[i, j])*(1 - self.Spec[j]))

      r = np.append(p, 1-sum(p))
      prod = np.array(prod)
      # log-probability, gradients[3], probability
      return [np.log(prod.dot(r)), prod[0]-prod[3], prod[1]-prod[3], prod[2]-prod[3], prod.dot(r)]

  # ###Fisher information for a test pattern
  def getPatternFisherInfo(self, p, t):
      x = self.txset.copy().drop(columns=['Test.Type'])
      # Compute the Fisher information matrix, weighted
      I = np.diag([0.0, 0.0, 0.0])
      # if t in [1, 3, 5, 7]: #consider only cases with RT-PCR tests
      for i in range(3):
          for j in range(i,3):
              sx = x.loc[self.txset['Test.Type'] == t+1]
              v2 = 0
              for k in range(sx.shape[0]):
                  y = sx.iloc[k,].values.astype(str).tolist()
                  Q = self.qxp(y, p) #Q[2:4] are gradients. Q[5] is probability
                  v2 = v2 + ((Q[1 + i]) * Q[1 + j]) / (Q[4])
              I[i, j] = v2
              I[j, i] = v2
      return I

  ###Weighted Fisher information
  def getWeightedFisherInfo(self, wsub, p):
      #Initialise to uniform test patterns. Ensure interior point, to begin.
      ttype  = self.txset['Test.Type'].unique().tolist()
      w = np.append(0, np.array(wsub))
      I = np.diag([0.0, 0.0, 0.0])
      for t in range(1, len(ttype)):
          widx = t - 1
          I = I + (w[t] / self.patternCost[widx]) * self.getPatternFisherInfo(p, t)
      return I

max_stars_repo_name: cni-iisc/optimal-sero-survey-design

max_issues_repo_name: cni-iisc/optimal-sero-survey-design

max_forks_repo_name: cni-iisc/optimal-sero-survey-design