gauss.c

#include <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>

#ifndef COEFF_ARR_SIZE
#define COEFF_ARR_SIZE 8
#endif

float a0, a1, a2, a3, b1, b2, left_corner, right_corner;

float *gauss_coefficients(float sigma) {
  if (sigma < 0.5) {
    sigma = 0.5;
  }

  float a = exp(0.726 * 0.726) / sigma;
  float g1 = exp(-a);
  float g2 = exp(-2 * a);
  float k = (1 - g1) * (1 - g1) / (1 + 2 * a * g1 - g2);

  a0 = k;
  a1 = k * (a - 1) * g1;
  a2 = k * (a + 1) * g1;
  a3 = -k * g2;
  b1 = 2 * g1;
  b2 = -g2;
  left_corner = (a0 + a1) / (1 - b1 - b2);
  right_corner = (a2 + a3) / (1 - b1 - b2);

  float *coefficients = (float*) malloc(COEFF_ARR_SIZE * sizeof(float));

  coefficients[0] = a0;
  coefficients[1] = a1;
  coefficients[2] = a2;
  coefficients[3] = a3;
  coefficients[4] = b1;
  coefficients[5] = b2;
  coefficients[6] = left_corner;
  coefficients[7] = right_corner;

  return coefficients;
}

void convolve_mono_16(uint16_t* src, uint16_t* out, float* line,
		      float* coeff, size_t width, size_t height) {
  float coeff_a0, coeff_a1, coeff_b1, coeff_b2;
  
  uint16_t curr_src, prev_src, curr_out, prev_out, prev_prev_out;
  
  size_t src_index, out_index, line_index;

  for (size_t i = 0; i < height; ++i) {
    src_index = i * width;
    out_index = i;
    line_index = 0;

    prev_src = src[src_index];
    prev_prev_out = prev_src * coeff[6];
    prev_out = prev_prev_out;
    
    coeff_a0 = coeff[0];
    coeff_a1 = coeff[1];
    coeff_b1 = coeff[4];
    coeff_b2 = coeff[5];

    for (size_t j = 0; j < width; ++j) {
      curr_src = src[src_index];

      curr_out = curr_src * coeff_a0 +
	prev_src * coeff_a1 +
	prev_out * coeff_b1 +
	prev_prev_out * coeff_b2;

      prev_prev_out = prev_out;
      prev_out = curr_out;
      prev_src = curr_src;

      line[line_index] = prev_out;
      line_index++;
      src_index++;
    }

    src_index--;
    line_index--;
    out_index += height * (width - 1);

    // right to left
    prev_src = src[src_index];
    prev_prev_out = prev_src * coeff[7];
    prev_out = prev_prev_out;
    curr_src = prev_src;

    coeff_a0 = coeff[2];
    coeff_a1 = coeff[3];

    for (int64_t j = width - 1; j >= 0; --j) {
      curr_out = curr_src * coeff_a0 +
	prev_src * coeff_a1 +
	prev_out * coeff_b1 +
	prev_prev_out * coeff_b2;

      prev_prev_out = prev_out;
      prev_out = curr_out;

      prev_src = curr_src;
      curr_src = src[src_index];

      out[out_index] = line[line_index] + prev_out;

      src_index--;
      line_index--;
      out_index -= height;
    }
  }
}

int blur_mono_16(uint16_t *src, int src_length, int width,
		 int height, float radius) {
  if (!radius) return -1;

  uint16_t *out = (uint16_t*) malloc(src_length * sizeof(uint16_t));
  if (out == NULL) {
    return -1;
  }

  float *tmp_line = (float*) malloc(fmax(width, height) * sizeof(float));
  if (tmp_line == NULL) {
    return -1;
  }

  float *coeff = gauss_coefficients(radius);
  if (coeff == NULL) {
    return -1;
  }

  convolve_mono_16(src, out, tmp_line, coeff, width, height);
  convolve_mono_16(out, src, tmp_line, coeff, height, width);
  
  free(coeff);
  free(tmp_line);
  free(out);

  return 0;
}

int main() {
  size_t width = pow(2, 15);
  size_t height = pow(2, 15);
  float radius = 0.5;

  size_t src_size = width * height;

  uint16_t* src = (uint16_t*) malloc(src_size * sizeof(uint16_t));
  if (src == NULL) {
    return -1;
  }
  for (size_t i = 0; i < src_size; ++i) {
    src[i] = i;
  }

  int res = blur_mono_16(src, src_size, width, height, radius);

  free(src);
  return res;
}