'**SUMMARY**
Find szeilski or Harris features in the image.
Harris features correspond to Harris corner detection in the image.
Read more:
Harris Features: http://en.wikipedia.org/wiki/Corner_detection
szeliski Features: http://research.microsoft.com/en-us/um/people/szeliski/publications.htm
**PARAMETERS**
* *method* - Features type
* *threshold* - threshold val
**RETURNS**
A list of Feature objects corrseponding to the feature points.
**EXAMPLE**
>>> img = Image("corner_sample.png")
>>> fpoints = img.findFeatures("harris", 2000)
>>> for f in fpoints:
... f.draw()
>>> img.show()
**SEE ALSO**
:py:meth:`drawKeypointMatches`
:py:meth:`findKeypoints`
:py:meth:`findKeypointMatch`'

def findFeatures(self, method='szeliski', threshold=1000):

    try:
        import cv2
    except ImportError:
        logger.warning('OpenCV    >=    2.3.0    required')
        return None
    img = self.getGrayNumpyCv2()
    blur = cv2.GaussianBlur(img, (3, 3), 0)
    Ix = cv2.Sobel(blur, cv2.CV_32F, 1, 0)
    Iy = cv2.Sobel(blur, cv2.CV_32F, 0, 1)
    Ix_Ix = np.multiply(Ix, Ix)
    Iy_Iy = np.multiply(Iy, Iy)
    Ix_Iy = np.multiply(Ix, Iy)
    Ix_Ix_blur = cv2.GaussianBlur(Ix_Ix, (5, 5), 0)
    Iy_Iy_blur = cv2.GaussianBlur(Iy_Iy, (5, 5), 0)
    Ix_Iy_blur = cv2.GaussianBlur(Ix_Iy, (5, 5), 0)
    harris_thresh = (threshold * 5000)
    alpha = 0.06
    detA = ((Ix_Ix_blur * Iy_Iy_blur) - (Ix_Iy_blur ** 2))
    traceA = (Ix_Ix_blur + Iy_Iy_blur)
    feature_list = []
    if (method == 'szeliski'):
        harmonic_mean = (detA / traceA)
        for (j, i) in np.argwhere((harmonic_mean > threshold)):
            feature_list.append(Feature(self, i, j, ((i, j), (i, j), (i, j), (i, j))))
    elif (method == 'harris'):
        harris_function = (detA - ((alpha * traceA) * traceA))
        for (j, i) in np.argwhere((harris_function > harris_thresh)):
            feature_list.append(Feature(self, i, j, ((i, j), (i, j), (i, j), (i, j))))
    else:
        logger.warning('Invalid    method.')
        return None
    return feature_list
