/**
 * \file dirgrad.c
 * \brief This file contains the implementation of the Dir Grad Buffer creation, deletion and manipulation
 *
 *
 *  \ingroup corefiles
 *  \date Apr 14, 2009
 *  \author Alexander Pimenov
 *
 *  For the internal format for the DirGrad (or simply DG) buffer see the #PACK_DIRGRAD(X,Y) macro
 *
 *
 */

#include <stdlib.h>
#include "global.h"

#include "dirgrad.h"
#include "g12buffer.h"
#include "algebra.h"

#ifdef ASSERT
int dirgradBufferCounter = 0;
#endif

/**
 * \brief This function creates a dirgrad 12 bit buffer for the given width and height.
 *
 * In fact in the PC implementation a 16bit buffer is used. The created image is filled with 0 values
 * that represent the 0 magnitude and the 0 angle direction
 *
 * \attention You should not forget to destroy the created buffer with the call to dirGradBufferDelete(toDelete)
 *
 * \param h The height of the image to be created.
 * \param w The width of the image to be created.
 *
 * \return a newly created buffer.
 */
DirGradBuffer *dirGradBufferCreate(int h, int w)
{
	DirGradBuffer *toReturn = (DirGradBuffer *) dv_malloc (sizeof (DirGradBuffer));
	toReturn->h = h;
	toReturn->w = w;
	toReturn->data = (uint16_t *) dv_calloc(h * w , sizeof(uint16_t));

#ifdef ASSERT
	dirgradBufferCounter++;
#endif

	return toReturn;
}

/**
 * \brief This function deletes a DirGrad scale 12 bit buffer.
 *
 * \param[in] toDelete a buffer to be freed
 */
void dirGradBufferDelete(DirGradBuffer *toDelete)
{
	if (toDelete->data != NULL)
		dv_free(toDelete->data);
	dv_free(toDelete);
#ifdef ASSERT
	dirgradBufferCounter++;
#endif

}

/**
 * \brief This function creates a new DirGrad function form the given 12 bit greyscale buffer.
 *
 * This function computes the magnitude and the direction of the gradient and packs them into 12 bits.
 *
 * Sobel derivatives filters are used
 *
 *
 * \f$ P=\pmatrix{
 *       -1 & 0  & 1\cr
 *       -2 & 0  & 2\cr
 *       -1 & 0  & 1}
 * \f$
 * 		and
 * \f$ P=\pmatrix{
 *       -1 & 2 & -1\cr
 *        0 & 0 &  0\cr
 *        1 & 2 &  1}
 * \f$
 *
 *    The magnitude is then estimated as \f$\sqrt(Gx^2 + Gy^2)\f$ and the direction is \f$atan2(Gx, Gy)\f$
 *
 * \attention You should not forget to destroy the created buffer with the call to dirGradBufferDelete(toDelete)
 *
 */

DirGradBuffer *dirGradBufferCreateFromG12(G12Buffer *src)
{
	int i,j;
	int h = src->h;
	int w = src->w;
	DirGradBuffer* toReturn = dirGradBufferCreate(h,w);


	for (i = 1; i < h - 1; i++)
	{
		for (j = 1; j < w -1; j++ )
		{
			// TODO: Optimize this for optimal accessing

			// DX and DY will range in [ -2^12 ; 2^12 ]
			int16_t dx =    ( - src->data[(i - 1) * w + j - 1] + src->data[(i - 1) * w + j + 1]
			          + 2 * ( - src->data[      i * w + j - 1] + src->data[      i * w + j + 1])
			                  - src->data[(i + 1) * w + j - 1] + src->data[(i + 1) * w + j + 1]) / 4;

			int16_t dy =    ( - src->data[(i - 1) * w + j - 1] + src->data[(i + 1) * w + j - 1]
					  + 2 * ( - src->data[(i - 1) * w + j    ] + src->data[(i + 1) * w + j    ])
							  - src->data[(i - 1) * w + j + 1] + src->data[(i + 1) * w + j + 1]) / 4;
			// 7 bits for magnitude (12 - 8) * 2 - 1 = 7
			int magnitude = getMagnitude7bit(dx,dy);
			// 5 bits for gradient
			int direction = getDirection5bit(dx,dy);
			toReturn->data[i * w + j] = PACK_DIRGRAD(direction, magnitude);
		}
	}

	return toReturn;
}
