""" Tensorflow implementation of the face detection / alignment algorithm found at
https://github.com/kpzhang93/MTCNN_face_detection_alignment
"""
# MIT License
# 
# Copyright (c) 2016 David Sandberg
# 
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# 
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# 
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from six import string_types, iteritems
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
#from math import floor
import cv2
import os

def layer(op):
    '''Decorator for composable network layers.'''

    def layer_decorated(self, *args, **kwargs):
        # Automatically set a name if not provided.
        name = kwargs.setdefault('name', self.get_unique_name(op.__name__))
        # Figure out the layer inputs.
        if len(self.terminals) == 0:
            raise RuntimeError('No input variables found for layer %s.' % name)
        elif len(self.terminals) == 1:
            layer_input = self.terminals[0]
        else:
            layer_input = list(self.terminals)

        # Setup the layer and return it 
        layer_output = op(self, layer_input, *args, **kwargs)

        # Add to layer LUT.
        self.layers[name] = layer_output
        # This output is now the input for the next layer.
        self.feed(layer_output)
        # Return self for chained calls.

        # print("Added layer is: {0}".format(name))
        # print("layer_input {0}".format(layer_input))
        # print("arguments")
        # for arg in args:
            # print(arg)
        # print("kwargs: {0}".format(kwargs))
        # print(op)
        # print("layer_output: {0}", layer_output)
        # print("layer_output.get_shape(): {0}".format(layer_output.get_shape()))
        return self

    return layer_decorated

class Network(object):

    def __init__(self, inputs, trainable=True):
        # The input nodes for this network
        self.inputs = inputs
        # The current list of terminal nodes
        self.terminals = []
        # Mapping from layer names to layers
        self.layers = dict(inputs)
        # If true, the resulting variables are set as trainable
        self.trainable = trainable

        # print(self.layers) 
        # print(self.inputs)
        self.setup()

    def setup(self):
        '''Construct the network. '''
        raise NotImplementedError('Must be implemented by the subclass.')

    def load(self, data_path, session, ignore_missing=False):
        '''Load network weights.
        data_path: The path to the numpy-serialized network weights
        session: The current TensorFlow session
        ignore_missing: If true, serialized weights for missing layers are ignored.
        '''
        data_dict = np.load(data_path, encoding='latin1').item() #pylint: disable=no-member
        for op_name in data_dict:
            with tf.variable_scope(op_name, reuse=True):
                for param_name, data in iteritems(data_dict[op_name]):
                    try:
                        var = tf.get_variable(param_name)
                        session.run(var.assign(data))
                    except ValueError:
                        if not ignore_missing:
                            raise

    def feed(self, *args):
        '''Set the input(s) for the next operation by replacing the terminal nodes.
        The arguments can be either layer names or the actual layers.
        '''
        '''The purpose of this layer is to 
        '''
        assert len(args) != 0
        self.terminals = []
        for fed_layer in args:
            # print("fed_layer in feed is: {0}".format(fed_layer))
            # check if the argument is a string
            if isinstance(fed_layer, string_types):
                # if it is then check if it is in our layers dictionary
                try:
                    fed_layer = self.layers[fed_layer]
                except KeyError:
                    raise KeyError('Unknown layer name fed: %s' % fed_layer)
            # then append it to the terminal (which is the end of the network)
            self.terminals.append(fed_layer)
        return self

    def get_output(self):
        '''Returns the current network output.'''
        return self.terminals[-1]

    def get_unique_name(self, prefix):
        '''Returns an index-suffixed unique name for the given prefix.
        This is used for auto-generating layer names based on the type-prefix.
        '''
        ident = sum(t.startswith(prefix) for t, _ in self.layers.items()) + 1
        return '%s_%d' % (prefix, ident)

    def make_var(self, name, shape):
        '''Creates a new TensorFlow variable.'''
        return tf.get_variable(name, shape, trainable=self.trainable)

    def validate_padding(self, padding):
        '''Verifies that the padding is one of the supported ones.'''
        assert padding in ('SAME', 'VALID')

    # equivalent to conv = layer(conv)
    @layer
    def conv(self,
             inp,
             k_h,
             k_w,
             c_o,
             s_h,
             s_w,
             name,
             relu=True,
             padding='SAME',
             group=1,
             biased=True):
        # print("inp: {0}".format(inp))
        # print("k_h: {0}".format(k_h))
        # print("k_w: {0}".format(k_w))
        # print("c_o: {0}".format(c_o))
        # print("s_h: {0}".format(s_h))
        # print("s_w: {0}".format(s_w))
        # Verify that the padding is acceptable
        self.validate_padding(padding)
        # Get the number of channels in the input
        c_i = int(inp.get_shape()[-1])
        # Verify that the grouping parameter is valid
        assert c_i % group == 0
        assert c_o % group == 0
        # Convolution for a given input and kernel
        convolve = lambda i, k: tf.nn.conv2d(i, k, [1, s_h, s_w, 1], padding=padding)
        with tf.variable_scope(name) as scope:
            kernel = self.make_var('weights', shape=[k_h, k_w, c_i // group, c_o])
            # This is the common-case. Convolve the input without any further complications.
            output = convolve(inp, kernel)
            # Add the biases
            if biased:
                biases = self.make_var('biases', [c_o])
                output = tf.nn.bias_add(output, biases)
            if relu:
                # ReLU non-linearity
                output = tf.nn.relu(output, name=scope.name)
            return output

    @layer
    def prelu(self, inp, name):
        with tf.variable_scope(name):
            i = int(inp.get_shape()[-1])
            alpha = self.make_var('alpha', shape=(i,))

            # print("name in prelu is: {0}".format(name))
            # print("number of channels in input layer for prelu layer is: {0}".format(i))
            # print("name of alpha is: {0}".format(alpha.name))

            # Note here that the output just gets the negative values from the input vector,
            # multiplied by some trainable alpha value!
            output = tf.nn.relu(inp) + tf.multiply(alpha, -tf.nn.relu(-inp))
        return output

    # equivalent to max_pool = layer(max_pool)
    @layer
    def max_pool(self, inp, k_h, k_w, s_h, s_w, name, padding='SAME'):
        self.validate_padding(padding)
        return tf.nn.max_pool(inp,
                              ksize=[1, k_h, k_w, 1],
                              strides=[1, s_h, s_w, 1],
                              padding=padding,
                              name=name)

    @layer
    def fc(self, inp, num_out, name, relu=True):
        with tf.variable_scope(name):
            input_shape = inp.get_shape()
            if input_shape.ndims == 4:
                # The input is spatial. Vectorize it first.
                dim = 1
                for d in input_shape[1:].as_list():
                    dim *= int(d)
                feed_in = tf.reshape(inp, [-1, dim])
            else:
                feed_in, dim = (inp, input_shape[-1].value)
            weights = self.make_var('weights', shape=[dim, num_out])
            biases = self.make_var('biases', [num_out])
            op = tf.nn.relu_layer if relu else tf.nn.xw_plus_b
            fc = op(feed_in, weights, biases, name=name)
            return fc


    """
    Multi dimensional softmax,
    refer to https://github.com/tensorflow/tensorflow/issues/210
    compute softmax along the dimension of target
    the native softmax only supports batch_size x dimension
    """
    @layer
    def softmax(self, target, axis, name=None):
        max_axis = tf.reduce_max(target, axis, keep_dims=True)
        target_exp = tf.exp(target-max_axis)
        normalize = tf.reduce_sum(target_exp, axis, keep_dims=True)
        softmax = tf.div(target_exp, normalize, name)
        return softmax
    
class PNet(Network):
    def setup(self):
        (self.feed('data') # feeding the data into the network 
             .conv(3, 3, 10, 1, 1, padding='VALID', relu=False, name='conv1') #kernel_width, kernel_height, output_depth, stride_height, stride_width
             .prelu(name='PReLU1')
             .max_pool(2, 2, 2, 2, name='pool1')
             .conv(3, 3, 16, 1, 1, padding='VALID', relu=False, name='conv2')
             .prelu(name='PReLU2')
             .conv(3, 3, 32, 1, 1, padding='VALID', relu=False, name='conv3')
             .prelu(name='PReLU3')
             .conv(1, 1, 2, 1, 1, relu=False, name='conv4-1')
             .softmax(3,name='prob1'))
        (self.feed('PReLU3') #pylint: disable=no-value-for-parameter
             .conv(1, 1, 4, 1, 1, relu=False, name='conv4-2'))
        
class RNet(Network):
    def setup(self):
        (self.feed('data').conv(3, 3, 28, 1, 1, padding='VALID', relu=False, name='conv1')
             .prelu(name='prelu1')
             .max_pool(3, 3, 2, 2, name='pool1')
             .conv(3, 3, 48, 1, 1, padding='VALID', relu=False, name='conv2')
             .prelu(name='prelu2')
             .max_pool(3, 3, 2, 2, padding='VALID', name='pool2')
             .conv(2, 2, 64, 1, 1, padding='VALID', relu=False, name='conv3')
             .prelu(name='prelu3')
             .fc(128, relu=False, name='conv4')
             .prelu(name='prelu4')
             .fc(2, relu=False, name='conv5-1')
             .softmax(1,name='prob1'))

        (self.feed('prelu4') #pylint: disable=no-value-for-parameter
             .fc(4, relu=False, name='conv5-2'))

class ONet(Network):
    def setup(self):
        (self.feed('data') #pylint: disable=no-value-for-parameter, no-member
             .conv(3, 3, 32, 1, 1, padding='VALID', relu=False, name='conv1')
             .prelu(name='prelu1')
             .max_pool(3, 3, 2, 2, name='pool1')
             .conv(3, 3, 64, 1, 1, padding='VALID', relu=False, name='conv2')
             .prelu(name='prelu2')
             .max_pool(3, 3, 2, 2, padding='VALID', name='pool2')
             .conv(3, 3, 64, 1, 1, padding='VALID', relu=False, name='conv3')
             .prelu(name='prelu3')
             .max_pool(2, 2, 2, 2, name='pool3')
             .conv(2, 2, 128, 1, 1, padding='VALID', relu=False, name='conv4')
             .prelu(name='prelu4')
             .fc(256, relu=False, name='conv5')
             .prelu(name='prelu5')
             .fc(2, relu=False, name='conv6-1')
             .softmax(1, name='prob1'))

        (self.feed('prelu5') #pylint: disable=no-value-for-parameter
             .fc(4, relu=False, name='conv6-2'))

        (self.feed('prelu5') #pylint: disable=no-value-for-parameter
             .fc(10, relu=False, name='conv6-3'))

def create_mtcnn(sess, model_path):
    # If there isn't a model path
    if not model_path:
        # Then get the model path, note here that __file__ simply the file name
        model_path,_ = os.path.split(os.path.realpath(__file__)) 
        # print(model_path, _)


    with tf.variable_scope('pnet'):
        data = tf.placeholder(tf.float32, (None,None,None,3), 'input')
        pnet = PNet({'data':data})
        pnet.load(os.path.join(model_path, 'det1.npy'), sess)
    with tf.variable_scope('rnet'):
        data = tf.placeholder(tf.float32, (None,24,24,3), 'input')
        rnet = RNet({'data':data})
        rnet.load(os.path.join(model_path, 'det2.npy'), sess)
    with tf.variable_scope('onet'):
        data = tf.placeholder(tf.float32, (None,48,48,3), 'input')
        onet = ONet({'data':data})
        onet.load(os.path.join(model_path, 'det3.npy'), sess)
    
    # Setup the lambda function to be called to run the predict of the image 
    pnet_fun = lambda img : sess.run(('pnet/conv4-2/BiasAdd:0', 'pnet/prob1:0'), feed_dict={'pnet/input:0':img})
    rnet_fun = lambda img : sess.run(('rnet/conv5-2/conv5-2:0', 'rnet/prob1:0'), feed_dict={'rnet/input:0':img})
    onet_fun = lambda img : sess.run(('onet/conv6-2/conv6-2:0', 'onet/conv6-3/conv6-3:0', 'onet/prob1:0'), feed_dict={'onet/input:0':img})
    return pnet_fun, rnet_fun, onet_fun

def detect_face(img, minsize, pnet, rnet, onet, threshold, factor):
    # im: input image
    # minsize: minimum of faces' size
    # pnet, rnet, onet: caffemodel
    # threshold: threshold=[th1 th2 th3], th1-3 are three steps's threshold
    # fastresize: resize img from last scale (using in high-resolution images) if fastresize==true
    factor_count=0
    total_boxes=np.empty((0,9))
    points=[]
    h=img.shape[0]
    w=img.shape[1]
    minl=np.amin([h, w]) # get the minimum length
    m=12.0/minsize
    minl=minl*m # minimum_width_of_picture * 12/minimum_size_of_face
    # create scale pyramid
    scales=[]
    while minl>=12:
        scales += [m*np.power(factor, factor_count)]
        minl = minl*factor
        factor_count += 1
    # print("scales are: {}".format(scales))
    # first stage
    # print("shape of image is: {}".format(img.shape))
    for j in range(len(scales)):
        scale=scales[j]
        hs=int(np.ceil(h*scale))
        ws=int(np.ceil(w*scale))
        im_data = imresample(img, (hs, ws))
        im_data = (im_data-127.5)*0.0078125 # What is this constant for?
        img_x = np.expand_dims(im_data, 0)
        img_y = np.transpose(img_x, (0,2,1,3)) # switching rgb channels or maybe points for bbox
        out = pnet(img_y) # conv4-2/BiasAdd:0 has an output depth of 4 for the bounding boxes
                          # prob1:0 has an output depth of 2 to determine whether or not it is a human
        out0 = np.transpose(out[0], (0,2,1,3)) # this is a 4-tuple vector with the regression coefficients for the bounding box
        out1 = np.transpose(out[1], (0,2,1,3)) # this is a 2-tuple vector with each element being a probability on the convoluted cells
        # print("out is {}".format(out))
        # print("out0 is {}".format(out0))
        # print("out0.shape is {}".format(out0.shape))
        # print("out1 is {}".format(out1))
        # print("out1.shape is {}".format(out1.shape))      
        # # print("out0[0,:,:,:] is {}".format(out0[0,:,:,:]))
        # # print("out0[0,:,:,:].shape is {}".format(out0[0,:,:,:].shape))
        # # print("out1[0:,:,:,1] is {}".format(out1[0,:,:,1]))
        # # print("out1[0:,:,:,1].shape is {}".format(out0[0,:,:,:1].shape))

        # boxes is  = np.hstack([q1, q2, np.expand_dims(score,1))
        # _ is = out0 (the regression coefficients for the bounding box)
        # where q1 and q2 are the coordinates of the top left and bottom right corner of an image
        boxes, _ = generateBoundingBox(out1[0,:,:,1].copy(), out0[0,:,:,:].copy(), scale, threshold[0])
        # print("boxes are: {}".format(boxes.astype(np.int))) 
        # print("boxes shape is: {}".format(boxes.shape)) # (num_of_bbox, 5) 
        # inter-scale nms
        pick = nms(boxes.copy(), 0.5, 'Union') 

        # After each type of different scaling used on the image to find the bboxes, add the bboxes that are 
        # found to the total_boxes array
        if boxes.size>0 and pick.size>0:
            boxes = boxes[pick,:] # get the picked boxes from nms
            # total_boxes is a
            total_boxes = np.append(total_boxes, boxes, axis=0) # add them to the array
    numbox = total_boxes.shape[0]

    # print("Stage 1: Number of boxes after detecting bboxes in all scales are: {}".format(numbox))
    # print("total_boxes is: {}".format(total_boxes))
    # print("total_boxes.shape is: {}".format(total_boxes.shape))

    if numbox>0:
        # NMS again to refine the boxes
        pick = nms(total_boxes.copy(), 0.7, 'Union')
        total_boxes = total_boxes[pick,:]
        regw = total_boxes[:,2]-total_boxes[:,0] # width of all boxes
        regh = total_boxes[:,3]-total_boxes[:,1] # height of all boxes

        # find the corners of all the boxes
        qq1 = total_boxes[:,0]+total_boxes[:,5]*regw 
        qq2 = total_boxes[:,1]+total_boxes[:,6]*regh
        qq3 = total_boxes[:,2]+total_boxes[:,7]*regw
        qq4 = total_boxes[:,3]+total_boxes[:,8]*regh


        total_boxes = np.transpose(np.vstack([qq1, qq2, qq3, qq4, total_boxes[:,4]]))
        total_boxes = rerec(total_boxes.copy()) # change all the boxes into squares
        total_boxes[:,0:4] = np.fix(total_boxes[:,0:4]).astype(np.int32) # round down the coordinates for the corners
        dy, edy, dx, edx, y, ey, x, ex, tmpw, tmph = pad(total_boxes.copy(), w, h)

    numbox = total_boxes.shape[0]
    
    # print("Stage 1: after refinement of bboxes: {}".format(numbox))
    # print("total_boxes is: {}".format(total_boxes))
    # print("total_boxes.shape is: {}".format(total_boxes.shape))

    if numbox>0:
        # second stage
        tempimg = np.zeros((24,24,3,numbox)) # create a place to store the bounding boxes
        for k in range(0,numbox):
            tmp = np.zeros((int(tmph[k]),int(tmpw[k]),3)) # temporary empty image that is the size of a bbox
            tmp[dy[k]-1:edy[k],dx[k]-1:edx[k],:] = img[y[k]-1:ey[k],x[k]-1:ex[k],:] # set the tmp empty image's pixels 
            # plt.imshow(tmp)
            # image_save_file = os.path.join(os.getcwd(), "tmp.png") 
            # # print("Saving heat map image to: {}".format(image_save_file))
            # plt.savefig(image_save_file)
            # resize the image extracted by the bbox to 24, 2
            if tmp.shape[0]>0 and tmp.shape[1]>0 or tmp.shape[0]==0 and tmp.shape[1]==0:
                tempimg[:,:,:,k] = imresample(tmp, (24, 24))
            else:
                return np.empty()
        tempimg = (tempimg-127.5)*0.0078125 # normalize image
        tempimg1 = np.transpose(tempimg, (3,1,0,2)) # change to a format that is  (num of images, width of image, height of image, number of color channels)
        out = rnet(tempimg1) 
        out0 = np.transpose(out[0]) # vector of bounding boxes
        out1 = np.transpose(out[1]) # vector of probabilities
        # print("out0.shape is: {}".format(out0.shape)) # (4, num_of_bboxes) shape
        # print("out1.shape is: {}".format(out1.shape)) # (2, num_of_bboxes) shape
        # print("out0 is: {}".format(out0))
        # print("out1 is: {}".format(out1))
        score = out1[1,:] # get the probability that it is a face
        ipass = np.where(score>threshold[1]) # only get the indexes of the bboxes whereby the probability of said cell is over some threshold

        # Using the probabilities that a bbox is a face from rnet, filter out the bboxes that are of a lower probability to be a face
        total_boxes = np.hstack([total_boxes[ipass[0],0:4].copy(), np.expand_dims(score[ipass].copy(),1)])
        # print("total_boxes.shape after second stage is: {}".format(total_boxes.shape))
        # print("total_boxes after second stage is: {}".format(total_boxes))
        mv = out0[:,ipass[0]]
        if total_boxes.shape[0]>0:
            pick = nms(total_boxes, 0.7, 'Union')
            total_boxes = total_boxes[pick,:] # get the boxes
            total_boxes = bbreg(total_boxes.copy(), np.transpose(mv[:,pick])) # now using the regression to find boxes
            total_boxes = rerec(total_boxes.copy()) # turn boxes into squares
        # exit(0)

    numbox = total_boxes.shape[0]
    # print("Total number of bboxes found in second stage is: {}".format(numbox))
    # print("total_boxes is: {}".format(total_boxes))
    # print("total_boxes.shape is: {}".format(total_boxes.shape))
    if numbox>0:
        # third stage
        total_boxes = np.fix(total_boxes).astype(np.int32)
        dy, edy, dx, edx, y, ey, x, ex, tmpw, tmph = pad(total_boxes.copy(), w, h)
        tempimg = np.zeros((48,48,3,numbox))
        for k in range(0,numbox):
            tmp = np.zeros((int(tmph[k]),int(tmpw[k]),3))
            tmp[dy[k]-1:edy[k],dx[k]-1:edx[k],:] = img[y[k]-1:ey[k],x[k]-1:ex[k],:]
            if tmp.shape[0]>0 and tmp.shape[1]>0 or tmp.shape[0]==0 and tmp.shape[1]==0:
                tempimg[:,:,:,k] = imresample(tmp, (48, 48))
            else:
                return np.empty()
        tempimg = (tempimg-127.5)*0.0078125
        tempimg1 = np.transpose(tempimg, (3,1,0,2)) # change to a format that is  (num of images, width of image, height of image, number of color channels) 
        out = onet(tempimg1)

        out0 = np.transpose(out[0]) # regression coefficients for the bbox 
        out1 = np.transpose(out[1]) # facial feature points where the first 5 are the x coordinates and the last 5 are the corresponding y coordinates
        out2 = np.transpose(out[2]) # probability that bbox is a face
        # print("out0.shape is: {}".format(out0.shape)) # (4, num_of_bboxes) shape
        # print("out1.shape is: {}".format(out1.shape)) # (10, num_of_bboxes) shape
        # print("out2.shape is: {}".format(out2.shape)) # (2, num_of_bboxes) shape
        # print("out0 is: {}".format(out0))
        # print("out1 is: {}".format(out1))
        # print("out2 is: {}".format(out2))
        score = out2[1,:]
        points = out1
        ipass = np.where(score>threshold[2]) # find the indexes of the scores that is above the threshold probability
        points = points[:,ipass[0]] # get the facial points that are above the threshold probability 
        total_boxes = np.hstack([total_boxes[ipass[0],0:4].copy(), np.expand_dims(score[ipass].copy(),1)]) # get the bboxes

        # get the facial points that have a high probability to be face
        mv = out0[:,ipass[0]] 

        # calculate the width and height of all the bboxes
        w = total_boxes[:,2]-total_boxes[:,0]+1
        h = total_boxes[:,3]-total_boxes[:,1]+1

        # print("w is :{}".format(w))
        # print("np.tile(w, (5, 1)) is: {}".format(np.tile(w, (5, 1))))
        # print("points[0:5,:] is: {}".format(points[0:5,:]))
        # print("np.tile(total_boxes[:,0],(5, 1)) is: {}".format(np.tile(total_boxes[:,0], (5, 1))))

        # Scaling the points 
        # Note here that since points is a (10, num_of_bboxes) shape, each put needs to be
        # multiplied by the width so we create 5 copies of the width to multiply with points[0:5, :]
        # Then we add the value of the x-coordinate to the points 
        points[0:5,:] = np.tile(w,(5, 1))*points[0:5,:] + np.tile(total_boxes[:,0],(5, 1))-1
        points[5:10,:] = np.tile(h,(5, 1))*points[5:10,:] + np.tile(total_boxes[:,1],(5, 1))-1
        if total_boxes.shape[0]>0:
            total_boxes = bbreg(total_boxes.copy(), np.transpose(mv))
            pick = nms(total_boxes.copy(), 0.7, 'Min')
            total_boxes = total_boxes[pick,:]
            points = points[:,pick]


#        numbox = total_boxes.shape[0]
        # print("Total number of bboxes found in third stage is: {}".format(numbox))
        # print("total_boxes is: {}".format(total_boxes.astype(np.int)))
        # print("total_boxes.shape is: {}".format(total_boxes.shape))
#        draw_img = img.copy()
#        for bbox in total_boxes.astype(np.int):
#            x1,y1,x2,y2 = bbox[0], bbox[1], bbox[2], bbox[3]
#            cv2.rectangle(draw_img, (x1, y1), (x2, y2), (0,0,255), 5)
#        plt.imshow(draw_img)  
#        image_save_file = os.path.join(os.getcwd(), "bboxes.png")
#        plt.savefig(image_save_file)

    return total_boxes, points


def bulk_detect_face(images, detection_window_size_ratio, pnet, rnet, onet, threshold, factor):
    # im: input image
    # minsize: minimum of faces' size
    # pnet, rnet, onet: caffemodel
    # threshold: threshold=[th1 th2 th3], th1-3 are three steps's threshold [0-1]

    all_scales = [None] * len(images)
    images_with_boxes = [None] * len(images)

    for i in range(len(images)):
        images_with_boxes[i] = {'total_boxes': np.empty((0, 9))}

    # create scale pyramid
    for index, img in enumerate(images):
        all_scales[index] = []
        h = img.shape[0]
        w = img.shape[1]
        minsize = int(detection_window_size_ratio * np.minimum(w, h))
        factor_count = 0
        minl = np.amin([h, w])
        if minsize <= 12:
            minsize = 12

        m = 12.0 / minsize
        minl = minl * m
        while minl >= 12:
            all_scales[index].append(m * np.power(factor, factor_count))
            minl = minl * factor
            factor_count += 1

    # # # # # # # # # # # # #
    # first stage - fast proposal network (pnet) to obtain face candidates
    # # # # # # # # # # # # #

    images_obj_per_resolution = {}

    # TODO: use some type of rounding to number module 8 to increase probability that pyramid images will have the same resolution across input images

    for index, scales in enumerate(all_scales):
        h = images[index].shape[0]
        w = images[index].shape[1]

        for scale in scales:
            hs = int(np.ceil(h * scale))
            ws = int(np.ceil(w * scale))

            if (ws, hs) not in images_obj_per_resolution:
                images_obj_per_resolution[(ws, hs)] = []

            im_data = imresample(images[index], (hs, ws))
            im_data = (im_data - 127.5) * 0.0078125
            img_y = np.transpose(im_data, (1, 0, 2))  # caffe uses different dimensions ordering
            images_obj_per_resolution[(ws, hs)].append({'scale': scale, 'image': img_y, 'index': index})

    for resolution in images_obj_per_resolution:
        images_per_resolution = [i['image'] for i in images_obj_per_resolution[resolution]]
        outs = pnet(images_per_resolution)

        for index in range(len(outs[0])):
            scale = images_obj_per_resolution[resolution][index]['scale']
            image_index = images_obj_per_resolution[resolution][index]['index']
            out0 = np.transpose(outs[0][index], (1, 0, 2))
            out1 = np.transpose(outs[1][index], (1, 0, 2))

            boxes, _ = generateBoundingBox(out1[:, :, 1].copy(), out0[:, :, :].copy(), scale, threshold[0])

            # inter-scale nms
            pick = nms(boxes.copy(), 0.5, 'Union')
            if boxes.size > 0 and pick.size > 0:
                boxes = boxes[pick, :]
                images_with_boxes[image_index]['total_boxes'] = np.append(images_with_boxes[image_index]['total_boxes'],
                                                                          boxes,
                                                                          axis=0)

    for index, image_obj in enumerate(images_with_boxes):
        numbox = image_obj['total_boxes'].shape[0]
        if numbox > 0:
            h = images[index].shape[0]
            w = images[index].shape[1]
            pick = nms(image_obj['total_boxes'].copy(), 0.7, 'Union')
            image_obj['total_boxes'] = image_obj['total_boxes'][pick, :]
            regw = image_obj['total_boxes'][:, 2] - image_obj['total_boxes'][:, 0]
            regh = image_obj['total_boxes'][:, 3] - image_obj['total_boxes'][:, 1]
            qq1 = image_obj['total_boxes'][:, 0] + image_obj['total_boxes'][:, 5] * regw
            qq2 = image_obj['total_boxes'][:, 1] + image_obj['total_boxes'][:, 6] * regh
            qq3 = image_obj['total_boxes'][:, 2] + image_obj['total_boxes'][:, 7] * regw
            qq4 = image_obj['total_boxes'][:, 3] + image_obj['total_boxes'][:, 8] * regh
            image_obj['total_boxes'] = np.transpose(np.vstack([qq1, qq2, qq3, qq4, image_obj['total_boxes'][:, 4]]))
            image_obj['total_boxes'] = rerec(image_obj['total_boxes'].copy())
            image_obj['total_boxes'][:, 0:4] = np.fix(image_obj['total_boxes'][:, 0:4]).astype(np.int32)
            dy, edy, dx, edx, y, ey, x, ex, tmpw, tmph = pad(image_obj['total_boxes'].copy(), w, h)

            numbox = image_obj['total_boxes'].shape[0]
            tempimg = np.zeros((24, 24, 3, numbox))

            if numbox > 0:
                for k in range(0, numbox):
                    tmp = np.zeros((int(tmph[k]), int(tmpw[k]), 3))
                    tmp[dy[k] - 1:edy[k], dx[k] - 1:edx[k], :] = images[index][y[k] - 1:ey[k], x[k] - 1:ex[k], :]
                    if tmp.shape[0] > 0 and tmp.shape[1] > 0 or tmp.shape[0] == 0 and tmp.shape[1] == 0:
                        tempimg[:, :, :, k] = imresample(tmp, (24, 24))
                    else:
                        return np.empty()

                tempimg = (tempimg - 127.5) * 0.0078125
                image_obj['rnet_input'] = np.transpose(tempimg, (3, 1, 0, 2))

    # # # # # # # # # # # # #
    # second stage - refinement of face candidates with rnet
    # # # # # # # # # # # # #

    bulk_rnet_input = np.empty((0, 24, 24, 3))
    for index, image_obj in enumerate(images_with_boxes):
        if 'rnet_input' in image_obj:
            bulk_rnet_input = np.append(bulk_rnet_input, image_obj['rnet_input'], axis=0)

    out = rnet(bulk_rnet_input)
    out0 = np.transpose(out[0])
    out1 = np.transpose(out[1])
    score = out1[1, :]

    i = 0
    for index, image_obj in enumerate(images_with_boxes):
        if 'rnet_input' not in image_obj:
            continue

        rnet_input_count = image_obj['rnet_input'].shape[0]
        score_per_image = score[i:i + rnet_input_count]
        out0_per_image = out0[:, i:i + rnet_input_count]

        ipass = np.where(score_per_image > threshold[1])
        image_obj['total_boxes'] = np.hstack([image_obj['total_boxes'][ipass[0], 0:4].copy(),
                                              np.expand_dims(score_per_image[ipass].copy(), 1)])

        mv = out0_per_image[:, ipass[0]]

        if image_obj['total_boxes'].shape[0] > 0:
            h = images[index].shape[0]
            w = images[index].shape[1]
            pick = nms(image_obj['total_boxes'], 0.7, 'Union')
            image_obj['total_boxes'] = image_obj['total_boxes'][pick, :]
            image_obj['total_boxes'] = bbreg(image_obj['total_boxes'].copy(), np.transpose(mv[:, pick]))
            image_obj['total_boxes'] = rerec(image_obj['total_boxes'].copy())

            numbox = image_obj['total_boxes'].shape[0]

            if numbox > 0:
                tempimg = np.zeros((48, 48, 3, numbox))
                image_obj['total_boxes'] = np.fix(image_obj['total_boxes']).astype(np.int32)
                dy, edy, dx, edx, y, ey, x, ex, tmpw, tmph = pad(image_obj['total_boxes'].copy(), w, h)

                for k in range(0, numbox):
                    tmp = np.zeros((int(tmph[k]), int(tmpw[k]), 3))
                    tmp[dy[k] - 1:edy[k], dx[k] - 1:edx[k], :] = images[index][y[k] - 1:ey[k], x[k] - 1:ex[k], :]
                    if tmp.shape[0] > 0 and tmp.shape[1] > 0 or tmp.shape[0] == 0 and tmp.shape[1] == 0:
                        tempimg[:, :, :, k] = imresample(tmp, (48, 48))
                    else:
                        return np.empty()
                tempimg = (tempimg - 127.5) * 0.0078125
                image_obj['onet_input'] = np.transpose(tempimg, (3, 1, 0, 2))

        i += rnet_input_count

    # # # # # # # # # # # # #
    # third stage - further refinement and facial landmarks positions with onet
    # # # # # # # # # # # # #

    bulk_onet_input = np.empty((0, 48, 48, 3))
    for index, image_obj in enumerate(images_with_boxes):
        if 'onet_input' in image_obj:
            bulk_onet_input = np.append(bulk_onet_input, image_obj['onet_input'], axis=0)

    out = onet(bulk_onet_input)

    out0 = np.transpose(out[0])
    out1 = np.transpose(out[1])
    out2 = np.transpose(out[2])
    score = out2[1, :]
    points = out1

    i = 0
    ret = []
    for index, image_obj in enumerate(images_with_boxes):
        if 'onet_input' not in image_obj:
            ret.append(None)
            continue

        onet_input_count = image_obj['onet_input'].shape[0]

        out0_per_image = out0[:, i:i + onet_input_count]
        score_per_image = score[i:i + onet_input_count]
        points_per_image = points[:, i:i + onet_input_count]

        ipass = np.where(score_per_image > threshold[2])
        points_per_image = points_per_image[:, ipass[0]]

        image_obj['total_boxes'] = np.hstack([image_obj['total_boxes'][ipass[0], 0:4].copy(),
                                              np.expand_dims(score_per_image[ipass].copy(), 1)])
        mv = out0_per_image[:, ipass[0]]

        w = image_obj['total_boxes'][:, 2] - image_obj['total_boxes'][:, 0] + 1
        h = image_obj['total_boxes'][:, 3] - image_obj['total_boxes'][:, 1] + 1
        points_per_image[0:5, :] = np.tile(w, (5, 1)) * points_per_image[0:5, :] + np.tile(
            image_obj['total_boxes'][:, 0], (5, 1)) - 1
        points_per_image[5:10, :] = np.tile(h, (5, 1)) * points_per_image[5:10, :] + np.tile(
            image_obj['total_boxes'][:, 1], (5, 1)) - 1

        if image_obj['total_boxes'].shape[0] > 0:
            image_obj['total_boxes'] = bbreg(image_obj['total_boxes'].copy(), np.transpose(mv))
            pick = nms(image_obj['total_boxes'].copy(), 0.7, 'Min')
            image_obj['total_boxes'] = image_obj['total_boxes'][pick, :]
            points_per_image = points_per_image[:, pick]

            ret.append((image_obj['total_boxes'], points_per_image))
        else:
            ret.append(None)

        i += onet_input_count
    return ret


# function [boundingbox] = bbreg(boundingbox,reg)
def bbreg(boundingbox,reg):
    # calibrate bounding boxes
    if reg.shape[1]==1:
        reg = np.reshape(reg, (reg.shape[2], reg.shape[3]))

    w = boundingbox[:,2]-boundingbox[:,0]+1
    h = boundingbox[:,3]-boundingbox[:,1]+1
    b1 = boundingbox[:,0]+reg[:,0]*w
    b2 = boundingbox[:,1]+reg[:,1]*h
    b3 = boundingbox[:,2]+reg[:,2]*w
    b4 = boundingbox[:,3]+reg[:,3]*h
    boundingbox[:,0:4] = np.transpose(np.vstack([b1, b2, b3, b4 ]))
    return boundingbox
 
def generateBoundingBox(imap, reg, scale, t):
    # use heatmap to generate bounding boxes
    # imap is a heat map
    # reg is a matrix that has the x0, y0, width, height of a bounding box
    # t is some threshold probability
    # print("imap is {}".format(imap))
    # print("reg is {}".format(reg))
    stride=2
    cellsize=12

    imap = np.transpose(imap) # Why do we need to transpose the heat map?

    # Plot and save the heat map
#    plt.imshow(imap, cmap = 'hot')
#    image_save_file = os.path.join(os.getcwd(), "imap.png") 
    # # print("Saving heat map image to: {}".format(image_save_file))
#    plt.savefig(image_save_file)

    # print("reg is: {}".format(reg))
    # print("reg.shape is {}".format(np.array(reg).shape))

    # (dx1, dy1) is one corner
    # (dx2, dy2) is another corner
    dx1 = np.transpose(reg[:,:,0]) 
    dy1 = np.transpose(reg[:,:,1])
    dx2 = np.transpose(reg[:,:,2])
    dy2 = np.transpose(reg[:,:,3])
    # print("dx1.shape is: {}".format(dx1.shape))
    # print("dx1 is: {}".format(dx1))
    # print("dy1 is: {}".format(dy1))
    # print("dx2 is: {}".format(dx2))
    # print("dy2 is: {}".format(dy2))
    y, x = np.where(imap >= t) # return the indexes where imap >= t
    # print("y.shape is: {}".format(y.shape)) 
    if y.shape[0]==1:
        dx1 = np.flipud(dx1)
        dy1 = np.flipud(dy1)
        dx2 = np.flipud(dx2)
        dy2 = np.flipud(dy2)
    score = imap[(y,x)] # these are the probability values that above some threshold t for the convoluted cells
    # print("dx1[(y, x)].shape is: {}".format(dx1[(y, x)].shape))
    # print("dx1[(y, x)] is: {}".format(dx1[(y, x)]))
    # print("dy1[(y, x)] is: {}".format(dy1[(y, x)]))    
    # print("dx2[(y, x)] is: {}".format(dx2[(y, x)]))
    # print("dy2[(y, x)] is: {}".format(dy2[(y, x)]))

    # # print("score is: {}".format(score))
    reg = np.transpose(np.vstack([ dx1[(y,x)], dy1[(y,x)], dx2[(y,x)], dy2[(y,x)] ]))
    # print("new reg is: {}".format(reg)) # regression coefficients for bounding boxes that are past the threshold
    if reg.size==0:
        reg = np.empty((0,3))
    bb = np.transpose(np.vstack([y,x])) # these are the bounding boxes
    # print("bb is: {}".format(bb))
    # print("bb.shape is: {}".format(bb.shape))
    ## TODO: Understand how the image is scaled back to the original size to find the coordinates of the bboxes in 
    # the original image
    # print("scale is: {}".format(scale))
    # print("stride is: {}".format(stride))
    # print("cellsize is: {}".format(cellsize))
    q1 = np.fix((stride*bb+1)/scale) # rounds down any float to the nearest integer
    q2 = np.fix((stride*bb+cellsize-1+1)/scale) # rounds down any float to the nearest integer

    # q1 and q2 are the top left and bottom right corners coordinates in the real image
    # print("q1 is: {}".format(q1.astype(np.int32)))
    # print("q2 is: {}".format(q2.astype(np.int32)))
    # # print("np.expand_dims(score, 1) is: {}".format(np.expand_dims(score, 1)))
    # print("reg is: {}".format(reg))
    boundingbox = np.hstack([q1, q2, np.expand_dims(score,1), reg])
    # print("boundingbox.shape is : {}".format(boundingbox.shape))
    # print("boundingbox is: {}".format(boundingbox))
    return boundingbox, reg
 
# function pick = nms(boxes,threshold,type)
def nms(boxes, threshold, method):
    # boxes is the 9 dim input vector with the first 4 being the coordinates for the corners 
    # of the bounding boxes

    # threshold is the area that is union-ed
    if boxes.size==0:
        return np.empty((0,3))
    x1 = boxes[:,0] 
    y1 = boxes[:,1]
    x2 = boxes[:,2]
    y2 = boxes[:,3]
    s = boxes[:,4] # score of the bounding box
    # print("s is: {}".format(s))
    area = (x2-x1+1) * (y2-y1+1)
    # print("area is : {}".format(area))
    I = np.argsort(s)
    # print("I is: {}".format(I))
    pick = np.zeros_like(s, dtype=np.int16)
    # print("pick is: {}".format(pick))
    counter = 0
    # print(I.size)
    while I.size>0:
        i = I[-1] # pick the box that has the highest probability to be a face
        # print("I[-1] is: {}".format(i))
        pick[counter] = i # add the index of that box into the pick array
        counter += 1
        idx = I[0:-1] # get the rest of the boxes
        # print("idx is: {}".format(idx))

        # find the largest (x, y) coordinates for the start (top left) of the
        # bounding box and the smallest (x, y) coordinates for the end (bottom right)
        # for the end of the bounding box, the idea here is to pick the SMALLEST box

        # we are comparing the box that has the highest probability to be the face of a human
        # with all the other remaining boxes, and then getting the smallest box

        # idea here is that we want to get all the boxes surrounding the current box that is 
        # selected

        xx1 = np.maximum(x1[i], x1[idx]) 
        yy1 = np.maximum(y1[i], y1[idx])
        xx2 = np.minimum(x2[i], x2[idx])
        yy2 = np.minimum(y2[i], y2[idx])
    
        # Calculate the width and height for all the boxes
        w = np.maximum(0.0, xx2-xx1+1) # note here that xx2-xx1 will be negative if the boxes don't overlap
        h = np.maximum(0.0, yy2-yy1+1) # note here that yy2-yy1 will be negative if the boxes don't overlap

        # note here that inter stands for intersection
        inter = w * h # most of the values here will be 0 as only overlapping bounding boxes will have a value above 0
        # print("inter is: {}".format(inter))

        # Note here that o stands for overlap, and inter is the area of intersection
        if method is 'Min':
            o = inter / np.minimum(area[i], area[idx]) # this is the ratio of overlapping as inter
        else:
            o = inter / (area[i] + area[idx] - inter)
        # print("o is: {}".format(o))
        # print("I before thresholding: {}".format(I.shape))
        I = I[np.where(o<=threshold)] # now get keep the boxes that aren't overlapped as they might be new faces
        # print("I after thresholding: {}".format(I.shape))
    pick = pick[0:counter]
    return pick

# function [dy edy dx edx y ey x ex tmpw tmph] = pad(total_boxes,w,h)
def pad(total_boxes, w, h):
    # print("w is : {}".format(w))
    # print("h is : {}".format(h))
    # compute the padding coordinates (pad the bounding boxes to square)
    tmpw = (total_boxes[:,2]-total_boxes[:,0]+1).astype(np.int32) # width of all boxes
    tmph = (total_boxes[:,3]-total_boxes[:,1]+1).astype(np.int32) # height of all boxes
    numbox = total_boxes.shape[0]

    # values used to extract the pixels from the original image inside the bbox window
    dx = np.ones((numbox), dtype=np.int32) # this is used for padding
    dy = np.ones((numbox), dtype=np.int32) # this is used for padding 
    edx = tmpw.copy().astype(np.int32)
    edy = tmph.copy().astype(np.int32)

    # these are the actual values of the bbox in the original image
    x = total_boxes[:,0].copy().astype(np.int32)  
    y = total_boxes[:,1].copy().astype(np.int32)
    ex = total_boxes[:,2].copy().astype(np.int32)
    ey = total_boxes[:,3].copy().astype(np.int32)

    # print("x is: {}".format(x))
    # print("y is: {}".format(y))
    # print("ex is {}".format(ex)) 
    # print("ey is {}".format(ey))

    tmp = np.where(ex>w) # find the boxes where the x coordinate is greater than the picture
    # print("tmp is:{}".format(tmp))

    # updating the width so that the width is measured as: width of picture - x0,
    # instead of x1 - x0 (where x1 is outside the image dimension)
    edx.flat[tmp] = np.expand_dims(-ex[tmp]+w+tmpw[tmp],1)  
    # print(x[tmp], y[tmp], ex[tmp], ey[tmp])
    # print("ex[tmp] is: {}".format(ex[tmp]))
    # print("tmpw[tmp] is: {}".format(tmpw[tmp]))
    # print("edx.flat[tmp] is: {}".format(edx.flat[tmp]))
    # print("edx[tmp] is: {}".format(edx[tmp]))
    ex[tmp] = w  # setting the max value of the x coordinate of the corner to the max value of 
                 # the picture
    
    tmp = np.where(ey>h)
    # print("tmp is:{}".format(tmp))
    # updating the height so that the width is measured as: height of picture - y0,
    # instead of y1 - y0 (where y1 is outside the image dimension) 
    edy.flat[tmp] = np.expand_dims(-ey[tmp]+h+tmph[tmp],1) 
    # print("edy.flat[tmp] is: {}".format(edy.flat[tmp]))
    ey[tmp] = h 

    # finding the starting position of the bboxes which are less than 1
    tmp = np.where(x<1)
    # print("tmp is:{}".format(tmp))
    dx.flat[tmp] = np.expand_dims(2-x[tmp],1)
    # print("x[tmp] is: {}".format(x[tmp]))
    # print("dx.flat[tmp] is: {}".format(dx.flat[tmp]))
    x[tmp] = 1

    tmp = np.where(y<1)
    # print("tmp is:{}".format(tmp))
    dy.flat[tmp] = np.expand_dims(2-y[tmp],1) 
    # print("dy.flat[tmp] is: {}".format(dy.flat[tmp]))
    y[tmp] = 1
    return dy, edy, dx, edx, y, ey, x, ex, tmpw, tmph

# function [bboxA] = rerec(bboxA)
def rerec(bboxA):
    # convert bboxA to square
    h = bboxA[:,3]-bboxA[:,1]
    w = bboxA[:,2]-bboxA[:,0]
    l = np.maximum(w, h)
    # print("l is: {}".format(l))
    bboxA[:,0] = bboxA[:,0]+w*0.5-l*0.5
    bboxA[:,1] = bboxA[:,1]+h*0.5-l*0.5
    # print("bboxA[:, 2:4] is {}".format(bboxA[:, 2:4]))
    # print("np.transpose(np.tile(l, (2,1)) is: {}".format(np.transpose(np.tile(l,(2,1)))))
    bboxA[:,2:4] = bboxA[:,0:2] + np.transpose(np.tile(l,(2,1)))
    # print()
    return bboxA

def imresample(img, sz):
    im_data = cv2.resize(img, (sz[1], sz[0]), interpolation=cv2.INTER_AREA) #@UndefinedVariable
    return im_data

    # This method is kept for debugging purpose
#     h=img.shape[0]
#     w=img.shape[1]
#     hs, ws = sz
#     dx = float(w) / ws
#     dy = float(h) / hs
#     im_data = np.zeros((hs,ws,3))
#     for a1 in range(0,hs):
#         for a2 in range(0,ws):
#             for a3 in range(0,3):
#                 im_data[a1,a2,a3] = img[int(floor(a1*dy)),int(floor(a2*dx)),a3]
#     return im_data

