///////////////////////////////////////////////////////////////////////////////
// gradient.cpp
// ============
// 1st order derivative for horizontal and vertical direction.
// It computes horizontal and vertical gradient vectors of an image using 
// central difference.
// deltaX(x,y) = I(x+1, y) - I(x-1, y)
// deltaY(x,y) = I(x, y+1) - I(x, y-1)
//
// The gradients can be negative values, so the gradient output must be a 
// signed array. For example, for 8-bit char output, the range of values are
// from -128 to +127. If you want to store the results into "unsigned char", 
// use gradientNorm() instead. All values are shifted by 128 in gradientNorm().
// As a result, the normalized gradient image looks like with emboss effect.
//
// Mostly, the magnitude of the gragients are useful for image processing. It 
// is computed using;
// Mag(x,y) = sqrt(deltaX(x,y)*deltaX(x,y) + deltaY(x,y)*deltaY(x,y))
//
// Or, we can approximates by adding absolute values.
// Mag(x,y) ~ abs(deltaX(x,y)) + abs(deltaY(x,y))
//
//  AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
// CREATED: 2006-10-30
// UPDATED: 2006-10-30
///////////////////////////////////////////////////////////////////////////////

#include <cstring>
#include <cmath>
#include "gradient.h"

// function prototypes ////////////////////////////////////////////////////////
void gradientXY(const unsigned char* src, int width, int height, char* deltaX, char* deltaY);
void gradientXY(const unsigned char* src, int width, int height, int* deltaX, int* deltaY);

void gradientXYMag(const unsigned char* src, int width, int height, char* deltaX, char* deltaY, unsigned char* mag);
void gradientXYMag(const unsigned char* src, int width, int height, int* deltaX, int* deltaY, int* mag);

void gradientXYNorm(const unsigned char* src, int width, int height, unsigned char* deltaX, unsigned char* deltaY);
void gradientXYMagNorm(const unsigned char* src, int width, int height, unsigned char* deltaX, unsigned char* deltaY, unsigned char* mag);



///////////////////////////////////////////////////////////////////////////////
// compute gradient (first order derivative x and y) using central difference
// delta(t) = [f(t+1) - f(t-1)] / 2
// "char" version, the range is -128 ~ +127
///////////////////////////////////////////////////////////////////////////////
bool dip::gradient(const unsigned char* src, int width, int height, char* deltaX, char* deltaY, unsigned char* magnitude)
{
    if(!src || !deltaX || !deltaY)
        return false;

    // fill zeros at the first and last row
    memset(deltaX, 0, width);
    memset(deltaY, 0, width);
    memset(&deltaX[(height-1)*width], 0, width);
    memset(&deltaY[(height-1)*width], 0, width);
    // clear first and last column here
    int index = width - 1;
    for(int i = 1; i < height; ++i)             // loop for (height-1)
    {
        deltaX[index] = deltaX[index+1] = 0;    // the last column of current row and the first column of the next row
        deltaY[index] = deltaY[index+1] = 0;
        index += width;
    }

    // start to compute mag and grad
    if(magnitude)
    {
        // clear edges to zeros
        memset(magnitude, 0, width);
        memset(&magnitude[(height-1)*width], 0, width);
        index = width - 1;
        for(int i = 1; i < height; ++i)             // loop for (height-1)
        {
            magnitude[index] = magnitude[index+1] = 0;
            index += width;
        }

        // compute gradients and magnitude
        gradientXYMag(src, width, height, deltaX, deltaY, magnitude);
    }
    else
    {
        // compute gradients only
        gradientXY(src, width, height, deltaX, deltaY);
    }

    return true;
}



///////////////////////////////////////////////////////////////////////////////
// compute gradient (first order derivative x and y) using central difference
// delta(t) = f(t+1) - f(t-1)
// "int" version
///////////////////////////////////////////////////////////////////////////////
bool dip::gradient(const unsigned char* src, int width, int height, int* deltaX, int* deltaY, int* magnitude)
{
    if(!src || !deltaX || !deltaY)
        return false;

    // fill zeros at the first and last row
    memset(deltaX, 0, width);
    memset(deltaY, 0, width);
    memset(&deltaX[(height-1)*width], 0, width);
    memset(&deltaY[(height-1)*width], 0, width);
    // clear first and last column here
    int index = width - 1;
    for(int i = 1; i < height; ++i)             // loop for (height-1)
    {
        deltaX[index] = deltaX[index+1] = 0;    // the last column of current row and the first column of the next row
        deltaY[index] = deltaY[index+1] = 0;
        index += width;
    }

    // start to compute mag and grad
    if(magnitude)
    {
        // clear edges to zeros
        memset(magnitude, 0, width);
        memset(&magnitude[(height-1)*width], 0, width);
        index = width - 1;
        for(int i = 1; i < height; ++i)             // loop for (height-1)
        {
            magnitude[index] = magnitude[index+1] = 0;
            index += width;
        }

        // compute gradients and magnitude
        gradientXYMag(src, width, height, deltaX, deltaY, magnitude);
    }
    else
    {
        // compute gradients only
        gradientXY(src, width, height, deltaX, deltaY);
    }

    return true;
}



///////////////////////////////////////////////////////////////////////////////
// compute normalized gradient
// delta(t) = [f(t+1) - f(t-1)] / 2 + 128
// "char" version, the range is 0 ~ 255
///////////////////////////////////////////////////////////////////////////////
bool dip::gradientNorm(const unsigned char* src, int width, int height,
                       unsigned char* deltaX, unsigned char* deltaY, unsigned char* magnitude)
{
    if(!src || !deltaX || !deltaY)
        return false;

    // fill zeros at the first and last row
    memset(deltaX, 0, width);
    memset(deltaY, 0, width);
    memset(&deltaX[(height-1)*width], 0, width);
    memset(&deltaY[(height-1)*width], 0, width);
    // clear first and last column here
    int index = width - 1;
    for(int i = 1; i < height; ++i)             // loop for (height-1)
    {
        deltaX[index] = deltaX[index+1] = 0;    // the last column of current row and the first column of the next row
        deltaY[index] = deltaY[index+1] = 0;
        index += width;
    }

    // start to compute mag and grad
    if(magnitude)
    {
        // clear edges to zeros
        memset(magnitude, 0, width);
        memset(&magnitude[(height-1)*width], 0, width);
        index = width - 1;
        for(int i = 1; i < height; ++i)             // loop for (height-1)
        {
            magnitude[index] = magnitude[index+1] = 0;
            index += width;
        }

        // compute gradients and magnitude
        gradientXYMagNorm(src, width, height, deltaX, deltaY, magnitude);
    }
    else
    {
        // compute gradients only
        gradientXYNorm(src, width, height, deltaX, deltaY);
    }

    return true;
}



///////////////////////////////////////////////////////////////////////////////
// compute magnitude of gradient(deltaX & deltaY) per pixel
///////////////////////////////////////////////////////////////////////////////
bool dip::gradientMagnitude(const int* deltaX, const int* deltaY, int w, int h, int* mag)
{
    if(!deltaX || !deltaY || !mag)
        return false;

    int i, j, t;

#if(1) // accurate computation (SLOW)
    t = 0;
    for(i = 0; i < h; ++i)
        for(j = 0; j < w; ++j, ++t)
            mag[t] = (int)(sqrtf((float)deltaX[t]*deltaX[t] + (float)deltaY[t]*deltaY[t]) + 0.5f);

#else // fast estimation using abs()
    t = 0;
    for(i = 0; i < h; ++i)
        for(j = 0; j < w; ++j, ++t)
            mag[t] = abs(deltaX[t]) + abs(deltaY[t]);

#endif

    return true;
}



///////////////////////////////////////////////////////////////////////////////
// do both horizontal and vertical gradient
// char version, the range is -128 ~+127
///////////////////////////////////////////////////////////////////////////////
void gradientXY(const unsigned char* src, int w, int h, char* dx, char* dy)
{
    int i, j, t;
    int iMax = h - 1;
    int jMax = w - 1;

    t = w + 1;                  // the starting index at I[1, 1]
    for(i = 1; i < iMax; ++i)
    {
        for(j = 1; j < jMax; ++j)
        {
            // -128 -> 0, 0 ->128, 127 -> 255
            dx[t] = (char)((src[t + 1] - src[t - 1]) >> 1);
            dy[t] = (char)((src[t + w] - src[t - w]) >> 1);
            ++t;    // next column at the same row
        }
        ++t; ++t;   // next row
    }
}


///////////////////////////////////////////////////////////////////////////////
// do both horizontal and vertical gradient, plus magnitude
// char version, the range of gradient is -128 ~ +127
///////////////////////////////////////////////////////////////////////////////
void gradientXYMag(const unsigned char* src, int w, int h, char* dx, char* dy, unsigned char* mag)
{
    int i, j, t;
    int iMax = h - 1;
    int jMax = w - 1;

    t = w + 1;                  // the starting index at I[1, 1]
    for(i = 1; i < iMax; ++i)
    {
        for(j = 1; j < jMax; ++j)
        {
            // -128 -> 0, 0 ->128, 127 -> 255
            dx[t] = (char)((src[t + 1] - src[t - 1]) >> 1);
            dy[t] = (char)((src[t + w] - src[t - w]) >> 1);
            //mag[t] = (unsigned char)(sqrtf((float)(dx[t]*dx[t]) + (float)(dy[t]*dy[t]));
            mag[t] = (unsigned char)(abs(dx[t]) + abs(dy[t]));
            //NOTE: It is possible dx and dy both are -128. In this case,
            // the absolute sum dx and dy will be exceeded max value in char.
            ++t;    // next column at the same row
        }
        ++t; ++t;   // next row
    }
}



///////////////////////////////////////////////////////////////////////////////
// do both horizontal and vertical gradient
// int version, the range is -255 ~ +255
///////////////////////////////////////////////////////////////////////////////
void gradientXY(const unsigned char* src, int w, int h, int* dx, int* dy)
{
    int i, j, t;
    int iMax = h - 1;
    int jMax = w - 1;
    t = w + 1;                  // the starting index at I[1, 1]
    for(i = 1; i < iMax; ++i)
    {
        for(j=1, t=0; j < jMax; ++j)
        {
            // central difference, delta(t) = f(t+1) - f(t-1)
            dx[t] = src[t + 1] - src[t - 1];
            dy[t] = src[t + w] - src[t - w];
            ++t;    // next column at the same row
        }
        ++t; ++t;   // next row
    }
}



///////////////////////////////////////////////////////////////////////////////
// do both horizontal and vertical gradient, plus magnitude
// int version, the range of gradient is -255 ~ +255
///////////////////////////////////////////////////////////////////////////////
void gradientXYMag(const unsigned char* src, int w, int h, int* dx, int* dy, int* mag)
{
    int i, j, t;
    int iMax = h - 1;
    int jMax = w - 1;

    t = w + 1;                  // the starting index at I[1, 1]
    for(i = 1; i < iMax; ++i)
    {
        for(j = 1; j < jMax; ++j)
        {
            // central difference, delta(t) = f(t+1) - f(t-1)
            dx[t] = src[t + 1] - src[t - 1];
            dy[t] = src[t + w] - src[t - w];
            //mag[t] = (int)(sqrtf((float)(dx[t]*dx[t]) + (float)(dy[t]*dy[t])));
            mag[t] = abs(dx[t]) + abs(dy[t]);
            ++t;    // next column at the same row
        }
        ++t; ++t;   // next row
    }
}



///////////////////////////////////////////////////////////////////////////////
// do normalized horizontal and vertical gradient by shifting 128
// char version, the range is 0 ~ 255
///////////////////////////////////////////////////////////////////////////////
void gradientXYNorm(const unsigned char* src, int w, int h, unsigned char* dx, unsigned char* dy)
{
    int i, j, t;
    int iMax = h - 1;
    int jMax = w - 1;

    t = w + 1;                  // the starting index at I[1, 1]
    for(i = 1; i < iMax; ++i)
    {
        for(j = 1; j < jMax; ++j)
        {
            // -128 -> 0, 0 ->128, 127 -> 255
            dx[t] = (unsigned char)(((src[t + 1] - src[t - 1]) >> 1) + 128);
            dy[t] = (unsigned char)(((src[t + w] - src[t - w]) >> 1) + 128);
            ++t;    // next column at the same row
        }
        ++t; ++t;   // next row
    }
}



///////////////////////////////////////////////////////////////////////////////
// do normalized gradient by shifting 128, plus magnitude 
// char version, the range of gradient is 0 ~ 255
///////////////////////////////////////////////////////////////////////////////
void gradientXYMagNorm(const unsigned char* src, int w, int h, unsigned char* dx, unsigned char* dy, unsigned char* mag)
{
    int i, j, t;
    int iMax = h - 1;
    int jMax = w - 1;

    t = w + 1;                  // the starting index at I[1, 1]
    for(i = 1; i < iMax; ++i)
    {
        for(j = 1; j < jMax; ++j)
        {
            // -128 -> 0, 0 ->128, 127 -> 255
            dx[t] = (unsigned char)(((src[t + 1] - src[t - 1]) >> 1) + 128);
            dy[t] = (unsigned char)(((src[t + w] - src[t - w]) >> 1) + 128);
            mag[t] = (unsigned char)(abs(dx[t]-128) + abs(dy[t]-128));
            //NOTE: It is possible dx and dy both are -128. In this case,
            // the absolute sum dx and dy will be exceeded max value in char.
            ++t;    // next column at the same row
        }
        ++t; ++t;   // next row
    }
}



