standard library functions for Halide

std.cpp

#include "std_try.h"
#include <math.h>

using namespace Halide;

#define PI 3.14159
/*
  Interpolations
*/

// nearest neighbor interpolation
Expr NNinterpolation(Func f, Expr x, Expr y){
    Expr newX = min(x - cast<int>(x), cast<int>(x) + 1 -x);
    Expr newY = min(y - cast<int>(y), cast<int>(y) + 1 -y);
    return f(newX, newY);
}

//bilinear interpolation
Expr interpolate(Func f, Expr x, Expr y){
    Expr newX = cast<int>(x);
    Expr newY = cast<int>(y);

    Expr weight1 = newX+1 - x;
    Expr weight2 = newY+1 - y;
    Expr inter = (f(newX, newY)*weight1 + f(newX+1, newY)*(1-weight1))*weight2 +\
        (f(newX, newY+1)*weight1 + f(newX+1, newY+1)*(1-weight1))*(1-weight2);
    return inter;
}

//bicubic interpolation
Expr getValue(Expr p0, Expr p1, Expr p2, Expr p3, Expr coo){
    return p1 + 0.5f * coo * (p2 - p0 + coo*(2.0f*p0 - 5.0f*p1 + 4.0f*p2 - p3 + coo*(3.0f * (p1 - p2) + p3 - p0)));
}

Expr bicubic(Func f, Expr x, Expr y){
    Expr one = x - cast<int>(x);
    Expr two = y - cast<int>(y);
    Expr a1 = getValue(f(cast<int>(x)-1, cast<int>(y)-1), f(cast<int>(x), cast<int>(y)-1), \
                       f(cast<int>(x)+1, cast<int>(y)-1), f(cast<int>(x)+2, cast<int>(y)-1), one);
    Expr a2 = getValue(f(cast<int>(x)-1, cast<int>(y)), f(cast<int>(x), cast<int>(y)), \
                       f(cast<int>(x)+1, cast<int>(y)), f(cast<int>(x)+2, cast<int>(y)), one);
    Expr a3 = getValue(f(cast<int>(x)-1, cast<int>(y)+1), f(cast<int>(x), cast<int>(y)+1), \
                       f(cast<int>(x)+1, cast<int>(y)+1), f(cast<int>(x)+2, cast<int>(y)+1), one);
    Expr a4 = getValue(f(cast<int>(x)-1, cast<int>(y)+2), f(cast<int>(x), cast<int>(y)+2), \
                       f(cast<int>(x)+1, cast<int>(y)+2), f(cast<int>(x)+2, cast<int>(y)+2), one);
    return getValue(a1, a2, a3, a4, two);
}

/*
  Derivatives and Integrals
*/

//slope along a X axis
Func derivativeX(Func f){
    Func out;
    Var x,y;
    out(x,y) = f(x,y) - f(x-1,y);
    return out;
}

//slope along a Y axis
Func derivativeY(Func f){
    Func out;
    Var x,y;
    out(x,y) = f(x,y) - f(x, y-1);
    return out;
}

Func integrateX(Func f, Expr width){
    Func out;
    RVar rx(1, width-1);
    Var x,y;
    out(x,y) = f(0,y);
    out(rx,y) = out(rx-1,y) + f(rx,y);
    return out;
}

Func integrateY(Func f, Expr height){
    Func out;
    RVar ry(1, height-1);
    Var x,y;
    out(x,y) = f(x,0);
    out(x,ry) = out(x,ry-1) + f(x,ry);
    return out;
}

Func integrate(Func f, Expr width, Expr height){
    Func integX = integrateX(f, width);
    Func out = integrateY(integX, height);
    return out;
}

/*
  Color space conversions
*/

Func grayscale(Func f) {
    Var x,y,c;
    Func out;
    out(x,y,c) = 0.3f*f(x,y,0)+0.59f*f(x,y,1)+0.11f*f(x,y,2);
    return out;
}

Func rgbToYuv(Func f){
    Var x,y,c;
    Func out;
    Expr red = .299f*f(x,y,0) +  0.587f*f(x,y,1) + 0.114f*f(x,y,2);
    Expr green = -0.14713f*f(x,y,0) -0.28886f*f(x,y,1) + 0.436f*f(x,y,2);
    Expr blue = 0.615f*f(x,y,0) -0.51499f*f(x,y,1) -0.10001f*f(x,y,2);
    out(x,y) = (red, green, blue);
    return out;
}

Func yuvToRgb(Func f){
    Var x,y,c;
    Func out;
    Expr red = 1.0f*f(x,y,0) + 1.13983f*f(x,y,2);
    Expr green = 1.0f*f(x,y,0) - 0.39465f*f(x,y,1) - 0.58060f*f(x,y,2);
    Expr blue = 1.0f*f(x,y,0) + 2.03211f*f(x,y,1);
    out(x,y) = (red, green, blue);
    return out;
}

Func rgbToXyz(Func f){
    Var x,y;
    Func out;
    Expr r = f(x,y,0);
    Expr g = f(x,y,1);
    Expr b = f(x,y,2);
    float c1 = 0.04045f;
    float c2 = 1.0f/12.92f;
    float c3 = 0.055f;
    float c4 = 1.0f/(1.0f + c3);
    float c5 = 2.4f; //gamma
    r = select(r <= c1, r*c2, pow((r+c3)*c4, c5));
    g = select(g <= c1, g*c2, pow((g+c3)*c4, c5));
    b = select(b <= c1, b*c2, pow((b+c3)*c4, c5));
    Expr X = 0.4124f * r + 0.3576f * g + 0.1805f * b;
    Expr Y = 0.2126f * r + 0.7152f * g + 0.0722f * b;
    Expr Z = 0.0193f * r + 0.1192f * g + 0.9505f * b;
    out(x,y) = (X, Y, Z);
    return out;
}

Func xyzToRgb(Func f){
    Var x,y,z;
    Func out;
    Expr X = f(x,y,0);
    Expr Y = f(x,y,1);
    Expr Z = f(x,y,2);
    Expr r = 3.2406f * X - 1.5372f * Y - 0.4986 * Z;
    Expr g = -0.9689f * X + 1.8758f * Y + 0.0415f * Z;
    Expr b = 0.0557f * X - 0.2040f * Y + 1.0570f * Z;
    float c1 = 0.0031308f;
    float c2 = 12.92f;
    float c3 = 1.055f;
    float c4 = 0.055f;
    float c5 = 1.0f/2.4f;
    r = select(r <= c1, c2*r, c3*pow(r, c5) - c4);
    g = select(g <= c1, c2*g, c3*pow(g, c5) - c4);
    b = select(b <= c1, c2*b, c3*pow(b, c5) - c4);
    out(x,y) = (r,g,b);
    return out;
}

Func xyzToLab(Func f){
    Var x,y,z;
    Func out;
    Expr X = f(x,y,0);
    Expr Y = f(x,y,1);
    Expr Z = f(x,y,2);
    float Xn = 1.0f/(0.412453f + 0.357580f + 0.180423f);
    float Yn = 1.0f/(0.212671f + 0.715160f + 0.072169f);
    float Zn = 1.0f/(0.019334f + 0.119193f + 0.950227f);
    Expr one = select(cast<float>(Y*Yn) >= 0.008856f, cast<float>(pow(Y*Yn, 1.0f/3.0f)), cast<float>(7.787f*Y*Yn + 16.0f/116.0f));
    Expr two = select(cast<float>(X*Xn) >= 0.008856f, cast<float>(pow(X*Xn, 1.0f/3.0f)), cast<float>(7.787f*X*Xn + 16.0f/116.0f));
    Expr three = select(cast<float>(Z*Zn) >= 0.008856f, cast<float>(pow(Z*Zn, 1.0f/3.0f)), cast<float>(7.787f*Z*Zn + 16.0f/116.0f));
    Expr L = 1.16f * one - 0.16f;
    Expr a = 5.0f * (two - one);
    Expr b = 2.0f * (one - three);
    out(x,y) = (L, a, b);
    return out;
}

Func rgbToLab(Func f){
    return xyzToLab(rgbToXyz(f));
}

Func labToXyz(Func f){
    Func out;
    Var x,y;
    float s = 6.0f/29.0f;
    float Xn = 0.412453f + 0.357580f + 0.180423f;
    float Yn = 0.212671f + 0.715160f + 0.072169f;
    float Zn = 0.019334f + 0.119193f + 0.950227f;
    Expr L = f(x,y,0);
    Expr a = f(x,y,1);
    Expr b = f(x,y,2);
    Expr fy = (L + 0.16f)/1.16f;
    Expr fx = fy + a/5.0f;
    Expr fz = fy - b/2.0f;
    Expr Y = select(fy > s, Yn*fy*fy*fy, (fy - 16.0f/116.0f)*3*s*s*Yn);
    Expr X = select(fx > s, Xn*fx*fx*fx, (fx - 16.0f/116.0f)*3*s*s*Xn);
    Expr Z = select(fz > s, Zn*fz*fz*fz, (fz - 16.0f/116.0f)*3*s*s*Zn);
    out(x,y) = (X,Y,Z);
    return out;
}

Func labToRgb(Func f){
    return xyzToRgb(labToXyz(f));
}

Func rgbToHsv(Func f){
    Var x,y;
    Func out;
    float mult = 1.0f/6.0f;
    Expr r = f(x,y,0);
    Expr g = f(x,y,1);
    Expr b = f(x,y,2);
    Expr minV = min(r, min(g, b));
    Expr maxV = max(r, max(g, b));
    Expr v = maxV;
    Expr delta = cast<float>(maxV - minV);

    Expr s = select(delta == 0.0f, 0.0f, delta/maxV);
    Expr h = select(delta == 0.0f, 0.0f, select(r == maxV, 0.0f + (g - b) / delta,\
                                          select(g == maxV, 2.0f + (b - r) / delta,\
                                                 4.0f + (r - g) / delta)));
    h *= mult;
    h = select(h<0.0f, h+1, h);
    out(x,y) = (h, s, v);
    return out;
}

Func hsvToRgb(Func in){
    Var x,y;
    Func out;
    Expr h = 6*in(x,y,0);
    Expr s = in(x,y,1);
    Expr v = in(x,y,2);
    Expr i = select(cast<int>(h) == 6, 5, cast<int>(h));
    Expr f = h - i;
    Expr p = v * (1 - s);
    Expr q = v * (1 - s * f);
    Expr u = v * (1 - s * (1 - f));
    out(x,y) = select(s==0.0f, (v,v,v),\
                      select(i == 0, (v, u, p), \
                             select(i == 1, (q, v, p), \
                                    select(i == 2, (p, v, u),\
                                           select(i == 3, (p, q, v),\
                                                  select(i == 4, (u, p, v),\
                                                         (v, p, q)))))));
    return out;
}

/*
  Image scaling with bicubic interpolation
*/
Func resample(Func f, Expr factor){
    Var x,y;
    Func out;
    Expr newX = x/cast<double>(factor);
    Expr newY = y/cast<double>(factor);
    out(x,y) = bicubic(f, newX, newY);
    return out;
}

/*
  Gradient
*/
Func gradient(Func f){
    Var x,y;
    Func output;
    output(x,y) = -f(x-1,y-1) + f(x+1, y-1) - 2*f(x-1, y) + 2*f(x+1,y) - f(x-1, y+1) + f(x+1, y+1);
    return output;
}

/*
   Convolution with a given kernel
*/
Func convolution(Func f, Func kernel, Expr kernel_width, Expr kernel_height){
    Var x,y;
    Func convolved;
    RVar kx(0,kernel_width), ky(0,kernel_height);
    convolved(x,y) += kernel(kx,ky)*f(x+kx-(kernel_width/2),y+ky-(kernel_height/2));
    return convolved;
}

/*
  Some filters
*/
Func gaussianBlur(Func f, const float sigma){
    Var x,y;
    Func gaussian, kernelHor, kernelVer, hor, out;
    gaussian(x) = exp(-x*x/(2*sigma*sigma));
    float norm = 1.0/sqrt(2*PI*sigma*sigma);
    int radius = ceil(3*sigma);
    RVar i(-radius, 2*radius+1);
    kernelHor(x,y) = (norm*gaussian(x-radius))/sum(norm*gaussian(i));
    hor = convolution(f, kernelHor, 2*radius+1, 1);
    kernelVer(x,y) = (norm*gaussian(y-radius))/sum(norm*gaussian(i));
    out = convolution(hor, kernelVer, 1, 2*radius+1);
    return out;
}

Func bilateral(Func clamped, float sigmaDomain, float sigmaRange){
    RVar dx(cast<int>(-2*sigmaDomain), cast<int>(4*sigmaDomain+1)), dy(-cast<int>(2*sigmaDomain), cast<int>(4*sigmaDomain+1));
    Func g1, bil, euclidian,g2;
    Var i,j,x,y;
    g1(x,y) = (1.0/(2*sigmaDomain*sigmaDomain*PI))*exp(-(x*x+y*y)/float(2*sigmaDomain*sigmaDomain));
    euclidian(x,y,i,j) =sqrt(pow(clamped(x,y,0)-clamped(x+i, y+j,0),2)+ \
                         pow(clamped(x,y,1)-clamped(x+i, y+j,1),2)+   \
                         pow(clamped(x,y,2)-clamped(x+i, y+j,2),2));

    g2(x,y,i,j)= (1.0/sqrt(2*sigmaRange*sigmaRange*PI))*exp(-(euclidian(x,y,i,j)*euclidian(x,y,i,j))/float(2*sigmaRange*sigmaRange));
    bil(x,y) += clamped(x+dx, y+dy)*g1(dx, dy)*g2(x,y,dx,dy)/sum(g1(dx,dy)*g2(x,y,dx,dy));

    return bil;
}

/*
  Image histogram
*/
Func histogram(Func f, int buckets, Expr width, Expr height){
    Func hist;
    Var x;
    RVar rx(0, width), ry(0, height);
    hist(f(rx,ry))++;
    return hist;
}

/*
  Image rotation with bicubic interpolation
*/
Func rotate(Func f, Expr angle, Expr centerX, Expr centerY){
    Var x,y;
    Func out;
    Expr newX = (x-centerX)*cos(angle) - sin(angle)*(y-centerY) + centerX;
    Expr newY = sin(angle)*(x-centerX) + cos(angle)*(y-centerY) + centerY;
    out(x,y) = bicubic(f, newX, newY);
    return out;
}

std.h

#ifndef JOVANA_STD_TRY_H
#define JOVANA_STD_TRY_H

#include <Halide.h>
/*
  Color space conversions
*/
Halide::Func grayscale(Halide::Func f);
Halide::Func rgbToYuv(Halide::Func f);
Halide::Func yuvToRgb(Halide::Func f);
Halide::Func rgbToXyz(Halide::Func f);
Halide::Func xyzToRgb(Halide::Func f);
Halide::Func xyzToLab(Halide::Func f);
Halide::Func labToXyz(Halide::Func f);
Halide::Func labToRgb(Halide::Func f);
Halide::Func rgbToHsv(Halide::Func f);
Halide::Func hsvToRgb(Halide::Func f);
Halide::Func rgbToLab(Halide::Func f);

Halide::Func integrateX(Halide::Func input, Halide::Expr width);
Halide::Func integrateY(Halide::Func input, Halide::Expr height);
Halide::Func derivativeX(Halide::Func f);
Halide::Func derivativeY(Halide::Func f);
Halide::Func integrate(Halide::Func f, Halide::Expr width, Halide::Expr height);

Halide::Func rotate(Halide::Func input, Halide::Expr angle, Halide::Expr centerX, Halide::Expr centerY);
Halide::Func resample(Halide::Func f, Halide::Expr factor);

Halide::Func histogram(Halide::Func, int buckets, Halide::Expr width, Halide::Expr height);

Halide::Func convolution(Halide::Func f, Halide::Func kernel, Halide::Expr width, Halide::Expr height);
Halide::Func gradient(Halide::Func);
Halide::Func gaussianBlur(Halide::Func f, const float sigma);
Halide::Func bilateral(Halide::Func f, float spatialSigma,  float rangeSigma);


#endif