2018-12-16 05:11:41 +08:00
|
|
|
// SPDX-License-Identifier: MPL-2.0
|
|
|
|
// Copyright (c) Yuxuan Shui <yshuiv7@gmail.com>
|
|
|
|
|
2019-01-01 19:35:59 +08:00
|
|
|
#include <assert.h>
|
2018-12-16 05:11:41 +08:00
|
|
|
#include <math.h>
|
|
|
|
|
2019-01-21 00:53:39 +08:00
|
|
|
#include "compiler.h"
|
2018-12-16 05:11:41 +08:00
|
|
|
#include "kernel.h"
|
2018-12-20 04:50:02 +08:00
|
|
|
#include "utils.h"
|
2018-12-16 05:11:41 +08:00
|
|
|
|
2019-01-01 08:15:51 +08:00
|
|
|
/// Sum a region convolution kernel. Region is defined by a width x height rectangle whose
|
|
|
|
/// top left corner is at (x, y)
|
|
|
|
double sum_kernel(const conv *map, int x, int y, int width, int height) {
|
|
|
|
double ret = 0;
|
2018-12-16 05:11:41 +08:00
|
|
|
|
2019-02-18 05:54:35 +08:00
|
|
|
// Compute sum of values which are "in range"
|
|
|
|
int xstart = x, xend = width + x;
|
|
|
|
if (xstart < 0) {
|
2019-01-01 08:15:51 +08:00
|
|
|
xstart = 0;
|
2019-02-18 05:54:35 +08:00
|
|
|
}
|
|
|
|
if (xend > map->w) {
|
|
|
|
xend = map->w;
|
|
|
|
}
|
|
|
|
int ystart = y, yend = height + y;
|
|
|
|
if (ystart < 0) {
|
2019-01-01 08:15:51 +08:00
|
|
|
ystart = 0;
|
2019-02-18 05:54:35 +08:00
|
|
|
}
|
|
|
|
if (yend > map->h) {
|
|
|
|
yend = map->h;
|
|
|
|
}
|
|
|
|
assert(yend >= ystart && xend >= xstart);
|
2019-01-01 19:35:59 +08:00
|
|
|
|
2019-02-18 05:54:35 +08:00
|
|
|
int d = map->w;
|
2019-01-01 19:35:59 +08:00
|
|
|
if (map->rsum) {
|
2019-02-18 05:54:35 +08:00
|
|
|
// See sum_kernel_preprocess
|
2019-01-01 19:35:59 +08:00
|
|
|
double v1 = xstart ? map->rsum[(yend - 1) * d + xstart - 1] : 0;
|
|
|
|
double v2 = ystart ? map->rsum[(ystart - 1) * d + xend - 1] : 0;
|
2019-02-18 05:54:35 +08:00
|
|
|
double v3 = (xstart && ystart) ? map->rsum[(ystart - 1) * d + xstart - 1] : 0;
|
2019-01-01 19:35:59 +08:00
|
|
|
return map->rsum[(yend - 1) * d + xend - 1] - v1 - v2 + v3;
|
|
|
|
}
|
|
|
|
|
2019-01-01 08:15:51 +08:00
|
|
|
for (int yi = ystart; yi < yend; yi++) {
|
|
|
|
for (int xi = xstart; xi < xend; xi++) {
|
2019-01-01 19:35:59 +08:00
|
|
|
ret += map->data[yi * d + xi];
|
2018-12-16 05:11:41 +08:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2019-01-01 08:15:51 +08:00
|
|
|
return ret;
|
|
|
|
}
|
2018-12-16 05:11:41 +08:00
|
|
|
|
2019-01-01 08:15:51 +08:00
|
|
|
double sum_kernel_normalized(const conv *map, int x, int y, int width, int height) {
|
|
|
|
double ret = sum_kernel(map, x, y, width, height);
|
|
|
|
if (ret < 0) {
|
|
|
|
ret = 0;
|
|
|
|
}
|
|
|
|
if (ret > 1) {
|
|
|
|
ret = 1;
|
|
|
|
}
|
|
|
|
return ret;
|
2018-12-16 05:11:41 +08:00
|
|
|
}
|
|
|
|
|
2018-12-20 07:41:58 +08:00
|
|
|
static double attr_const gaussian(double r, double x, double y) {
|
2018-12-16 05:11:41 +08:00
|
|
|
// Formula can be found here:
|
|
|
|
// https://en.wikipedia.org/wiki/Gaussian_blur#Mathematics
|
|
|
|
// Except a special case for r == 0 to produce sharp shadows
|
|
|
|
if (r == 0)
|
|
|
|
return 1;
|
|
|
|
return exp(-0.5 * (x * x + y * y) / (r * r)) / (2 * M_PI * r * r);
|
|
|
|
}
|
|
|
|
|
|
|
|
conv *gaussian_kernel(double r) {
|
|
|
|
conv *c;
|
|
|
|
int size = r * 2 + 1;
|
|
|
|
int center = size / 2;
|
|
|
|
double t;
|
|
|
|
|
|
|
|
c = cvalloc(sizeof(conv) + size * size * sizeof(double));
|
2019-02-18 05:54:35 +08:00
|
|
|
c->w = c->h = size;
|
2019-01-01 19:35:59 +08:00
|
|
|
c->rsum = NULL;
|
2018-12-16 05:11:41 +08:00
|
|
|
t = 0.0;
|
|
|
|
|
|
|
|
for (int y = 0; y < size; y++) {
|
|
|
|
for (int x = 0; x < size; x++) {
|
|
|
|
double g = gaussian(r, x - center, y - center);
|
|
|
|
t += g;
|
|
|
|
c->data[y * size + x] = g;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
for (int y = 0; y < size; y++) {
|
|
|
|
for (int x = 0; x < size; x++) {
|
|
|
|
c->data[y * size + x] /= t;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
return c;
|
|
|
|
}
|
|
|
|
|
2019-01-01 08:15:51 +08:00
|
|
|
/// preprocess kernels to make shadow generation faster
|
|
|
|
/// shadow_sum[x*d+y] is the sum of the kernel from (0, 0) to (x, y), inclusive
|
2019-02-18 05:54:35 +08:00
|
|
|
void sum_kernel_preprocess(conv *map) {
|
|
|
|
if (map->rsum) {
|
2019-01-01 19:35:59 +08:00
|
|
|
free(map->rsum);
|
2019-02-18 05:54:35 +08:00
|
|
|
}
|
2019-01-01 08:15:51 +08:00
|
|
|
|
2019-02-18 05:54:35 +08:00
|
|
|
auto sum = map->rsum = ccalloc(map->w * map->h, double);
|
2019-01-01 08:15:51 +08:00
|
|
|
sum[0] = map->data[0];
|
|
|
|
|
2019-02-18 05:54:35 +08:00
|
|
|
for (int x = 1; x < map->w; x++) {
|
2019-01-01 08:15:51 +08:00
|
|
|
sum[x] = sum[x - 1] + map->data[x];
|
|
|
|
}
|
|
|
|
|
2019-02-18 05:54:35 +08:00
|
|
|
const int d = map->w;
|
|
|
|
for (int y = 1; y < map->h; y++) {
|
2019-01-01 08:15:51 +08:00
|
|
|
sum[y * d] = sum[(y - 1) * d] + map->data[y * d];
|
2019-02-18 05:54:35 +08:00
|
|
|
for (int x = 1; x < map->w; x++) {
|
2019-01-01 08:15:51 +08:00
|
|
|
double tmp = sum[(y - 1) * d + x] + sum[y * d + x - 1] -
|
|
|
|
sum[(y - 1) * d + x - 1];
|
|
|
|
sum[y * d + x] = tmp + map->data[y * d + x];
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2018-12-16 05:11:41 +08:00
|
|
|
// vim: set noet sw=8 ts=8 :
|