diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py
index a9adc76b..e577ab19 100644
--- a/skimage/feature/texture.py
+++ b/skimage/feature/texture.py
@@ -13,16 +13,15 @@ from ._texture import (_glcm_loop,
 
 def greycomatrix(image, distances, angles, levels=None, symmetric=False,
                  normed=False):
-    """Calculate the grey-level co-occurrence matrix.
+    """ *** YOYO *** Calculate the grey-level co-occurrence matrix.
 
     A grey level co-occurrence matrix is a histogram of co-occurring
     greyscale values at a given offset over an image.
 
     Parameters
     ----------
-    image : array_like of uint8
-        Integer typed input image. The image will be cast to uint8, so
-        the maximum value must be less than 256.
+    image : array_like of uint8 or uint16.
+        Integer typed input image. The image will be cast to uint16.
     distances : array_like
         List of pixel pair distance offsets.
     angles : array_like
@@ -30,10 +29,8 @@ def greycomatrix(image, distances, angles, levels=None, symmetric=False,
     levels : int, optional
         The input image should contain integers in [0, levels-1],
         where levels indicate the number of grey-levels counted
-        (typically 256 for an 8-bit image). Be aware that the co-occurrence
-        matrix for every angle and every distance is of size levels x levels. 
-        Choosing a too large level might result in exceedingly large matrix.
-        If you have 16 bit or 32 bit images, consider binning. 
+        (typically 256 for an 8-bit image). This argument is required for
+        16-bit images and is typically the maximum of the image. 
     symmetric : bool, optional
         If True, the output matrix `P[:, :, d, theta]` is symmetric. This
         is accomplished by ignoring the order of value pairs, so both
@@ -52,7 +49,8 @@ def greycomatrix(image, distances, angles, levels=None, symmetric=False,
         `P[i,j,d,theta]` is the number of times that grey-level `j`
         occurs at a distance `d` and at an angle `theta` from
         grey-level `i`. If `normed` is `False`, the output is of
-        type uint32, otherwise it is float64.
+        type uint32, otherwise it is float64. The dimensions are:
+        levels x levels x number of distances x number of angles.
 
     References
     ----------
@@ -103,15 +101,25 @@ def greycomatrix(image, distances, angles, levels=None, symmetric=False,
     assert image.min() >= 0
 
     image_max = image.max()
-    if levels is None:
-        # if levels is not given, we assume that there are [0, image_max -1] levels to 
-        # be considered. There would be only zeros for all other rows and columns. 
-        levels = image_max + 1
+
+    #if image.dtype in [np.float, np.float16, np.float32, np.float]
+    if image.dtype == np.float:
+        raise ValueError("Float images are not supported by greycomatrix."
+                         "The image needs to be casted to uint8 or uint16.")
+    
+    # for 16 bit images, levels must be set.
+    if image.dtype != np.uint8:
+        if levels is None:
+            raise ValueError("For images other than 8 bit images, argument levels must be set to a user defined value."
+                             "A large value will produce a large output matrix, potentially causing memory problems."
+                             "In this case, one might consider binning the image prior to calling greycomatrix.")
         
+    if levels is None:    
+        levels = 256
+
     assert image_max < levels
 
-    # we cast to uint16 (because of fixed typing in cython) 
-    # this has no impact on the size of the co-occurrence matrix. 
+    # image is cast to uint16 (because of fixed typing in cython) 
     image = image.astype(np.uint16)
 
     distances = np.ascontiguousarray(distances, dtype=np.float64)