def testRigidTransformEstimation(inImg, level, dTheta, displacement, thr):
    left=ndimage.rotate(inImg, dTheta)
    right=ndimage.rotate(inImg, -dTheta)
    left=ndimage.affine_transform(left , np.eye(2), offset=-1*displacement)
    right=ndimage.affine_transform(right, np.eye(2), offset=displacement)
    
    rightPyramid=[i for i in transform.pyramid_gaussian(right, level)]
    leftPyramid=[i for i in transform.pyramid_gaussian(left, level)]
    sel=level
    beta=estimateRigidTransformation(leftPyramid[sel], rightPyramid[sel], 2.0*dTheta, thr)
    return beta

def testRigidTransformationMultiscale(dTheta, displacement, level):
    inImg=misc.imread('T2sample.png')[...,0]
    left=ndimage.rotate(inImg, -0.5*dTheta)#Rotate symmetricaly to ensure both are still the same size
    right=ndimage.rotate(inImg, 0.5*dTheta)#Rotate symmetricaly to ensure both are still the same size
    right=ndimage.affine_transform(right, np.eye(2), offset=-1*displacement)
    rightPyramid=[i for i in transform.pyramid_gaussian(right, level)]
    leftPyramid=[i for i in transform.pyramid_gaussian(left, level)]
    rcommon.plotPyramids(leftPyramid, rightPyramid)
    beta=estimateRigidTransformationMultiscale(leftPyramid, rightPyramid)
    print 180.0*beta[0]/np.pi, beta[1:3]
    return beta

def main():
    # construct the argument parser and parse the arguments
    ap = argparse.ArgumentParser()
    ap.add_argument("-i", "--image", required=True, help="Path to the image")
    ap.add_argument("-s", "--scale", type=float, default=1.5, help="scale factor size")
    args = vars(ap.parse_args())

    # load the image
    image = cv2.imread(args["image"])

    # METHOD #1: No smooth, just scaling.
    # loop over the image pyramid
    for (i, resized) in enumerate(pyramid(image, scale=args["scale"])):
        # show the resized image
        cv2.imshow("Layer {}".format(i + 1), resized)
        cv2.waitKey(0)

    # close all windows
    cv2.destroyAllWindows()

    # METHOD #2: Resizing + Gaussian smoothing.
    for (i, resized) in enumerate(pyramid_gaussian(image, downscale=2)):
        # if the image is too small, break from the loop
        if resized.shape[0] < 30 or resized.shape[1] < 30:
            break

        # show the resized image
        cv2.imshow("Layer {}".format(i + 1), resized)
        cv2.waitKey(0)

def save_pyramid () :
    global temp_line
    global pyramids
    global patchNum
    global total_patch
    global total_pyramid

    org_img = Image.open("%s/%s.jpg" %(base_path, temp_line), 'r' )
    
    org_img_name = "%s " %(temp_line)        # original image name
    # describ_file.write ( "%s\n" %org_img_name )  # original image name
    
    pyramids = list( pyramid_gaussian(org_img, downscale=math.sqrt(2) ) )
    for i in range(len(pyramids) ):
        if min( pyramids[i].shape[0], pyramids[i].shape[1] ) < 30 :
            del pyramids[i:]
            break
    
    for i in range( len (pyramids) ) :
        row = pyramids[i].shape[0]
        col = pyramids[i].shape[1]
        im_matrix = np.zeros([row, col, 3]).astype('uint8')
    
        for k in range(row):
            for j in range(col):
                im_matrix[k,j] = pyramids[i][k,j] * 255
    
        new_img = Image.fromarray(im_matrix)
          # new_img.save("%s/pyramid-%s.jpg" %(patch_path, i+total_patch) )
        new_img.save("%s/pyramid-%s.jpg" %(patch_path, i+total_pyramid) )
        # new_img.show()
    
        patchNum[i] = (row-30+1) * (col-30+1)                  # the number of patches
    total_pyramid = total_pyramid + len(pyramids)
    total_patch = total_patch + sum(patchNum)

    def gaussian_pyramid(self, n_levels=3, downscale=2, sigma=None, order=1, mode="reflect", cval=0):
        r"""
        Return the gaussian pyramid of this image. The first image of the
        pyramid will be the original, unmodified, image.

        Parameters
        ----------
        n_levels : int
            Number of levels in the pyramid. When set to -1 the maximum
            number of levels will be build.

            Default: 3

        downscale : float, optional
            Downscale factor.

            Default: 2

        sigma : float, optional
            Sigma for gaussian filter. Default is `2 * downscale / 6.0` which
            corresponds to a filter mask twice the size of the scale factor
            that covers more than 99% of the gaussian distribution.

            Default: None

        order : int, optional
            Order of splines used in interpolation of downsampling. See
            `scipy.ndimage.map_coordinates` for detail.

            Default: 1

        mode :  {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
            The mode parameter determines how the array borders are handled,
            where cval is the value when mode is equal to 'constant'.

            Default: 'reflect'

        cval : float, optional
            Value to fill past edges of input if mode is 'constant'.

            Default: 0

        Returns
        -------
        image_pyramid:
            Generator yielding pyramid layers as menpo image objects.
        """
        max_layer = n_levels - 1
        pyramid = pyramid_gaussian(
            self.pixels, max_layer=max_layer, downscale=downscale, sigma=sigma, order=order, mode=mode, cval=cval
        )

        for j, image_data in enumerate(pyramid):
            image = self.__class__(image_data)

            # rescale and reassign existent landmark
            image.landmarks = self.landmarks
            transform = UniformScale(downscale ** j, self.n_dims)
            transform.pseudoinverse.apply_inplace(image.landmarks)
            yield image

def test_build_gaussian_pyramid():
    rows, cols, dim = image.shape
    pyramid = pyramid_gaussian(image, downscale=2)

    for layer, out in enumerate(pyramid):
        layer_shape = (rows / 2 ** layer, cols / 2 ** layer, dim)
        assert_array_equal(out.shape, layer_shape)

  def get_face(self, do_rot=True, do_scale=True):
    frame = self.ig.getFrame()
    frame_pyramid = list(tf.pyramid_gaussian(frame, max_layer=NUM_PYR, downscale=2))

    scale_ssds = {}
    for i, face_pyramid in enumerate(self.scaled_face_pyramids):
      if not do_scale and i != 1:
          continue
      res = self.determine_best_shift(face_pyramid, frame_pyramid)
      best_i, best_j, best_ssd = res
      scale_ssds[i] = (1.0 / (best_ssd * self.scaled_weights[i]), best_i, best_j, np.array(face_pyramid[0].shape))
    if len(scale_ssds) == 3 or not do_scale:
      best_i, best_j = scale_ssds[1][1], scale_ssds[1][2]
    else:
      best_i, best_j = scale_ssds[0][1], scale_ssds[0][2]
    total = sum([v[0] for v in scale_ssds.values()])
    interp_shape = sum([v[0] / total * v[3] for v in scale_ssds.values()])

    rot_ssds = {}
    for i, face_pyramid in enumerate(self.rotated_face_pyramids):
      if not do_rot and i != 1:
          continue
      res = self.determine_best_shift(face_pyramid, frame_pyramid)
      rot_best_i, rot_best_j, best_ssd = res
      rot_ssds[i] = (1.0 / best_ssd, rot_best_i, rot_best_j, np.array(face_pyramid[0].shape))
    total = sum([v[0] for v in rot_ssds.values()])
    interp_rot = sum([v[0] / total * ROT_AMTS[k] for k, v in rot_ssds.items()])

    return best_i, best_j, frame, interp_shape, interp_rot

  def calibrate(self):
    self.ig = WebcamImageGetter()
    self.ig.start()
    self.init_interp_shape = None
    print "Place face 1 ft from camera. When face is visible, press Enter to continue."
    while True:
      frame = self.ig.getFrame()
      if frame is None:
        continue
      gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
      face_cascade = cv2.CascadeClassifier("haarcascades/haarcascade_frontalface_default.xml")
      faces = face_cascade.detectMultiScale(gray, scaleFactor=1.3, minNeighbors=5)
      for x, y, w, h in faces:
        cv2.rectangle(frame, (x, y), (x + w, y + h), color=(255, 0, 0), thickness=2)
      cv2.imshow("calibration", frame)
      if cv2.waitKey(1) & 0xFF == 10:
        cv2.destroyWindow("calibration")
        if len(faces) > 0:
          break
        else:
          print "No face detected."

    x, y, w, h = faces[0]
    num_pix = float(w*h)
    face_roi = frame[y:y+h, x:x+w]
    rotated_faces = [tf.rotate(face_roi, angle=rot_ang) for rot_ang in ROT_AMTS]
    self.rotated_face_pyramids = [list(tf.pyramid_gaussian(face, max_layer=NUM_PYR, downscale=2))
                                  for face in rotated_faces]
    scaled_faces = [tf.rescale(face_roi, scale=sc) for sc in RESCALING_FACTORS]
    self.scaled_face_pyramids = [list(tf.pyramid_gaussian(face, max_layer=NUM_PYR, downscale=2))
                                 for face in scaled_faces]
    # scaled_weights are used for scaled_faces
    self.scaled_weights = [num_pix / (sf.shape[0]*sf.shape[1]) for sf in scaled_faces]
    # we observed that the small detector is too strong, so we penalize it more
    self.scaled_weights[0] *= 1.5
    # w = f*Y/Z  -->  f = wZ/Y
    self.camera_f = w * START_FACE_DIST/AVERAGE_FACE_WIDTH
    self.start_center = np.array((x + w/2.0, y+h/2.0))
    self.w = w; self.h = h
    cv2.destroyWindow("calibration")
    cv2.waitKey(1)
    cv2.destroyWindow("calibration")
    cv2.waitKey(1)
    print "Tracking face...press Enter to quit."
    print "Red: close, green: far, blue: in between."

def gaussian_downsample(frames, pyramid_levels=4):
    nt = frames.shape[0]
    for ii, frame in enumerate(frames):
        pyr = transform.pyramid_gaussian(frame.astype(np.float))
        for jj in xrange(pyramid_levels + 1):
            ds = pyr.next()
        if ii == 0:
            out = np.empty((nt,) + ds.shape, dtype=np.float)
        out[ii] = ds
    return out

def getpyramidImage(image, d, fname):
	
	path = fname + "/t%03d"%d + "_GaussianPyramidLevel"

	for (i, resized) in enumerate(pyramid_gaussian(image, downscale=2)):
		#if resized.shape[0] < 30 or resized.shape[1] < 30:
		if i > 4 :
		 	break
		#imsave("t000_GaussianPyramidLevel%i.tif"%i,resized)
		vigra.impex.writeVolume(convert(resized),path +"%i.tif"%i,'')

def getGaussianPyramidOfList(imageList,amountOfLayers):
	listOfGaussiansPyramids = list()
	for i in range(1,amountOfLayers+1):
		currentLayer = list()
		for currentImage in imageList:
			gaussianImage = tuple(transform.pyramid_gaussian(currentImage,max_layer=i))
			currentLayer.append(gaussianImage[-1])
		listOfGaussiansPyramids.append(currentLayer)

	return listOfGaussiansPyramids
		

def create_image_pyramid(img, downscale_dim=2, pyramid_layer=3):
    # I put in some automatic values for the definition call. I need to check that they actually work.
    """ Create image pyramid"""


    pyramid = tuple(pyramid_gaussian(pyramid_in, downscale=downscale_dim))

    """ Check pyramid results """
    #Also need to put some test cases that check that downscale_dim and pyramid_layer are correct values.

    return(pyramid)

def augment_data(img):
  rotate_angles = [0, 45, 90, 135, 180, 225, 270, 315]
  scales = 4 # number of downsampling scales
  flip_flags = [True, False]
  cnt = 0
  output_imgs, output_filenames = {}, {}

  for f in flip_flags:
    if f:
      arr_img = util.PIL2array(img)
      # plt.imshow(arr_img)
      f_img = flip_image(arr_img)
      # plt.imshow(f_img)
      f_img = util.array2PIL(f_img)

      """
      # Optional: using affine transformation
      # shear by 180 degrees is equivalent to rotation by 180 degrees + flip.
      # So after that we rotate it another 180 degrees to get just the flip.
      shear = 180
      rotation = 180
      tform_augment = transform.AffineTransform(scale=(1, 1), rotation=np.deg2rad(rotation),
                                               shear=np.deg2rad(shear), translation=(0, 0))
      f_img = transform.warp(arr_img, tform_augment, mode='constant', order=3)
      plt.imshow(f_img)

      """
    else:
      f_img = img

    pyramid = tuple(transform.pyramid_gaussian(f_img, downscale=2))
    for p in xrange(scales):
      H, W, chs = pyramid[p].shape
      p_img = util.array2PIL(pyramid[p])
      #plt.imshow(pyramid[p])
      #p_img.show()
      for angle in rotate_angles:
        output = p_img.rotate(angle, expand=True)
        output = output.resize((58, 58))
        output_imgs[cnt] = output
        # output.show()
        """
        if f:
          output.save('samples/' + 'flipped'+ '_p' + str(p+1) + '_r' + str(angle) + '.jpg')
        else:
          output.save('samples/' + 'p' + str(p + 1) + '_r' + str(angle) + '.jpg')
        """
        if f:
          output_filenames[cnt] = 'flipped' + '_p' + str(p + 1) + '_r' + str(angle) + '.jpg'
        else:
          output_filenames[cnt] = 'p' + str(p + 1) + '_r' + str(angle) + '.jpg'

        cnt += 1
  return output_imgs, output_filenames

def testEstimateRotationMultiscale(dTheta, level):
    #inImg=misc.imread('stinkbug.png')[...,0]
    inImg=misc.imread('T2sample.png')[...,0]
    left=ndimage.rotate(inImg, -dTheta/2.0)
    right=ndimage.rotate(inImg, dTheta/2.0)
    rightPyramid=[i for i in transform.pyramid_gaussian(right, level)]
    leftPyramid=[i for i in transform.pyramid_gaussian(left, level)]
    angles=[]
    theta=estimateRotationMultiscale(leftPyramid, rightPyramid, 0, angles)
    angles=180*np.array(angles).reshape(len(angles),1)/np.pi
    xticks=[str(i) for i in range(level+1)[::-1]]
    plt.figure()
    plt.plot(angles)
    plt.xlabel('Scale')
    plt.ylabel('Estimated angle')
    plt.title('Global rotation [GT/Est.]: '+str(dTheta)+'/'+str(180*theta/np.pi)+' deg. Levels: '+str(level+1))
    plt.xticks(range(level+1), xticks)
    plt.grid()
    #printPyramids(leftPyramid, rightPyramid)
    print 'Estimated:\n', angles,' degrees'
    return 180*theta/np.pi

    def image_pyramid_down(image, downscale=1.5, min_size=(30, 30)):
        """
        Method to downscale images and yield scaled images down to the provided min_size
        :param image: Image to downscale
        :param downscale: Downscale factor
        :param min_size: Minimum image size
        :return: Generator with scaled images
        """
        for i, resized_image in enumerate(pyramid_gaussian(image, downscale=downscale)):
            if resized_image.shape[0] < min_size[1] or resized_image.shape[1] < min_size[0]:
                break

            yield i, resized_image

def test_collection_viewer():

    img = data.astronaut()
    img_collection = tuple(pyramid_gaussian(img))

    view = CollectionViewer(img_collection)
    make_key_event(48)

    view.update_index('', 2),
    assert_equal(view.image, img_collection[2])
    view.keyPressEvent(make_key_event(53))
    assert_equal(view.image, img_collection[5])
    view._format_coord(10, 10)

def kMeans(img):
	t0 = time.time()
	# apply kMeans, fit data, get histogram, get color bar
	org_img = img
	img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
	z = img.reshape((-1, 3))
	#print(z.shape)

	# image resize just for silhouetteCoeff
	# Crops images to 300x300, but loses accuracy
	# Try pyrDown (downsampling the images)
	"""ysize, xsize, chan = img.shape
	if ysize and xsize > 300:
		xoff = (xsize - 300) // 2
		yoff = (ysize - 300) // 2
		y = img[yoff:-yoff, xoff:-xoff]
	else:
		y = img
	y = y.reshape((-1, 3))
	print(y.shape)"""

	# downnsample images with gaussian smoothing
	if (img.shape[0] > 250 or img.shape[1] > 250):
		for (i, resized) in enumerate(pyramid_gaussian(org_img, downscale=2)):
			if resized.shape[0] < 100 or resized.shape[1] < 100:
				break
			org_img = resized
		#cv2.imshow("Layer {}".format(i + 1), resized)
		#print(org_img.shape)

	org_img = org_img.reshape((-1, 3))
	org_img = scale(org_img)
	#print(org_img.shape)
	# kmeans
	clt = KMeans(n_clusters = silhouetteCoeff(org_img), random_state = 42)
	clt.fit(z)
	#print(clt.cluster_centers_)

	hist = centroidHistogram(clt)
	bar = plotColors(hist, clt.cluster_centers_)
	print("Time including KMeans: ", time.time() - t0)
	#print("unique labels: ", np.unique(np.array(clt.labels_), axis=0))

	plt.figure(1)
	plt.axis("off")
	plt.subplot(211)
	plt.imshow(img)
	plt.subplot(212)
	plt.imshow(bar)
	plt.show()

def compute_gaussian_pyramid(img, min_size):
    h, w = img.shape[0:2]
    curr_size = np.min([h, w])
    levels = 0

    while curr_size > min_size:
        curr_size = np.floor(curr_size / 2.0)
        levels += 1

    img_pyr = list(pyramid_gaussian(img, max_layer=levels))

    img_pyr.reverse()  # smallest to largest

    assert np.min(img_pyr[1].shape[:2]) > min_size
    assert np.min(img_pyr[1].shape[:2]) <= np.min(img_pyr[-1].shape[:2])

    return img_pyr

def create_ns (tmp_imgpath, cnt_ns ) :
	global pyramids

	tmp_img = Image.open("%s/%s" %(coco_path, tmp_imgpath), 'r' )
	pyramids = list( pyramid_gaussian( tmp_img, downscale=math.sqrt(2) ) )

	for i in range ( len(pyramids) ):
		if min( pyramids[i].shape[0], pyramids[i].shape[1] ) < MinFace :
			del pyramids[i:]
			break

	# for j in range(4) :
	for j in range(36) :
		# creating random index
		img_index = random.randint(0, len(pyramids)-1 )
		tmp_patch_num = ( pyramids[img_index].shape[0] - 12 + 1) * ( pyramids[img_index].shape[1] - 12 + 1)
		rand_index = random.randint(0, tmp_patch_num)

		# x, y position decoding
		row_max = pyramids[img_index].shape[0]
		col_max = pyramids[img_index].shape[1]
		row = 0
		col = rand_index
		
		while ( col >= col_max - 12 +1 ) :
			row = row + 1
			col = col - (col_max-12+1)

		flag = 0
		# Rejecting Black and White image
		tmp_ns = pyramids[img_index][row:row+12, col:col+12]
		if not len(tmp_ns.shape)==3 :
			print " Gray Image. Skip "
			return 0

		# Rejecting Positive Samples
		scale_factor = math.sqrt(2)**img_index

		tmp_ns = pyramids[img_index][row:row+12, col:col+12]
		tmp_ns = Image.fromarray((tmp_ns*255.0).astype(np.uint8) )
		# tmp_ns = tmp_ns.resize( (12,12), Image.BICUBIC )
		tmp_ns = tmp_ns.resize( (12,12), Image.BILINEAR )
		tmp_ns.save("%s/ns-%s.jpg" %(ns_path, cnt_ns+j) )

	return 1

def scale_pyramid(im_path):

    detections = []
    downscale = 1.5    
    
    # The current scale of the image
    scale = 0

    im = imread(im_path, as_grey=True)
    	
    top_2 = []

    # Downscale the image and iterate
    for im_scaled in pyramid_gaussian(im, downscale=downscale):

        if im_scaled.shape[0] < min_wdw_sz[1] or im_scaled.shape[1] < min_wdw_sz[0]:
            break

        for (x, y, im_window) in sliding_window(im_scaled, min_wdw_sz, step_size):
            if im_window.shape[0] != min_wdw_sz[1] or im_window.shape[1] != min_wdw_sz[0]:
                continue

            # Calculate the HOG features
            fd = hog(im_window, orientations, pixels_per_cell, cells_per_block, visualize, normalize)
            pred = clf.predict(fd)
            dec_score = clf.decision_function(fd)

            if(dec_score[0][pred - 1] > 0.5):
            	new_tuple = (x, y, dec_score[0][pred - 1], int(min_wdw_sz[0]*(downscale**scale)),
                	int(min_wdw_sz[1]*(downscale**scale)), pred)

            	if len(top_2) < 2:
            		top_2.append(new_tuple)
            	else:
            		if new_tuple[2] > top_2[0][2] and pred != top_2[0][5]:
            			top_2[0] = new_tuple
            		elif new_tuple[2] > top_2[1][2] and pred != top_2[0][5]:
            			top_2[1] = new_tuple
            		else:
            			continue
            
        scale+=1.25

    return top_2