stromnov · March 3, 2021 13:01 · May 9, 2012
diff --git a/std.cpp b/std.cpp
@@ -0,0 +1,342 @@
+#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;
+}
diff --git a/std.h b/std.h
@@ -0,0 +1,37 @@
+#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