#include <stdlib.h>
#include <math.h>
#include "darray.h"
-#include "colorize.h"
+#include "multigrid.h"
+
+/* if defined, residual is calculated and loop is aborted, if no improvement
+ * is reached or the target resudual change is reached. */
+#define RESIDUAL
/* just process, but do not generate neighbor list and do not solve multigrid */
//#define TEST_NO_RENDER
/* test expanding of X, Y and Z */
//#define TEST_EXPAND_Z
-/* print error message b, if a is false */
-#define ASSERT(a, b) \
-{ if (!(a)) { fprintf(stderr, "%s %s() Error in line %d: %s\n", __FILE__, __func__, __LINE__, (b)); abort(); } }
-
-void ASSERT_FLOW(int w, int h, int k, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy)
+void ASSERT_FLOW(int w, int h, int k, const darray_t *flow, const darray_t *flow_i)
{
- if (flow_x) {
- ASSERT(darrayGetNumberOfDimensions(flow_x) == 3, "flow_x array must have 3 dimensions");
- ASSERT(w == darrayGetDimensions(flow_x)[0], "x-sizes of flow_x vectors are different");
- ASSERT(h == darrayGetDimensions(flow_x)[1], "y-sizes of flow_x vectors are different");
- ASSERT(k-1 == darrayGetDimensions(flow_x)[2], "z-sizes of flow_x vectors are not k-1");
- }
- if (flow_y) {
- ASSERT(darrayGetNumberOfDimensions(flow_y) == 3, "flow_y array must have 3 dimensions");
- ASSERT(w == darrayGetDimensions(flow_y)[0], "x-sizes of flow_y vectors are different");
- ASSERT(h == darrayGetDimensions(flow_y)[1], "y-sizes of flow_y vectors are different");
- ASSERT(k-1 == darrayGetDimensions(flow_y)[2], "z-sizes of flow_y vectors are not k-1");
+ if (flow) {
+ ASSERT(darrayGetNumberOfDimensions(flow) == 4, "flow array must have 4 dimensions");
+ ASSERT(w == darrayGetDimensions(flow)[0], "x-sizes of flow vectors are different");
+ ASSERT(h == darrayGetDimensions(flow)[1], "y-sizes of flow vectors are different");
+ ASSERT(k-1 == darrayGetDimensions(flow)[2], "z-sizes of flow vectors are not k-1");
+ ASSERT(2 == darrayGetDimensions(flow)[3], "number of flow cubes are not 2");
}
- if (flow_ix) {
- ASSERT(darrayGetNumberOfDimensions(flow_ix) == 3, "flow_ix array must have 3 dimensions");
- ASSERT(w == darrayGetDimensions(flow_ix)[0], "x-sizes of flow_ix vectors are different");
- ASSERT(h == darrayGetDimensions(flow_ix)[1], "y-sizes of flow_ix vectors are different");
- ASSERT(k-1 == darrayGetDimensions(flow_ix)[2], "z-sizes of flow_ix vectors are not k-1");
- }
- if (flow_iy) {
- ASSERT(darrayGetNumberOfDimensions(flow_iy) == 3, "flow_iy array must have 3 dimensions");
- ASSERT(w == darrayGetDimensions(flow_iy)[0], "x-sizes of flow_iy vectors are different");
- ASSERT(h == darrayGetDimensions(flow_iy)[1], "y-sizes of flow_iy vectors are different");
- ASSERT(k-1 == darrayGetDimensions(flow_iy)[2], "z-sizes of flow_iy vectors are not k-1");
+ if (flow_i) {
+ ASSERT(darrayGetNumberOfDimensions(flow_i) == 4, "flow_i array must have 4 dimensions");
+ ASSERT(w == darrayGetDimensions(flow_i)[0], "x-sizes of flow_i vectors are different");
+ ASSERT(h == darrayGetDimensions(flow_i)[1], "y-sizes of flow_i vectors are different");
+ ASSERT(k-1 == darrayGetDimensions(flow_i)[2], "z-sizes of flow_i vectors are not k-1");
+ ASSERT(2 == darrayGetDimensions(flow_i)[3], "number of flow_i cubes are not 2");
}
}
-#define min(a, b) ((a < b) ? a : b)
-#define max(a, b) ((a > b) ? a : b)
-
#define NEIGHBORS_2D 9
#define NEIGHBORS_3D 27
-struct neighbors_2d {
- int num; /* # of neighbors */
- int index[NEIGHBORS_2D]; /* index of neighbors in an array */
- double weight[NEIGHBORS_2D]; /* calculated weight of neighbors */
-};
-struct neighbors_3d {
- int num; /* # of neighbors */
- int index[NEIGHBORS_3D]; /* index of neighbors in an array */
- double weight[NEIGHBORS_3D]; /* calculated weight of neighbors */
-};
-
static struct neigh_off {
int x, y, z;
} neigh_off[NEIGHBORS_3D] = {
*
* src plane 4 -------> dst plane 2
*/
-static int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow_x, const darray_t *flow_y)
+int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow)
{
int x, y, z, xx, yy, w, h, ks, kd;
int i, pixles;
ASSERT(w == darrayGetDimensions(dst)[0], "x-size of source and destination must be equal");
ASSERT(h == darrayGetDimensions(dst)[1], "y-size of source and destination must be equal");
ASSERT(kd == ((ks+1)>>1), "z-size of destination must be half of the source (+an extra plane if source z-size is odd)");
- ASSERT_FLOW(w, h, ks, flow_x, flow_y, NULL, NULL);
+ ASSERT_FLOW(w, h, ks, flow, NULL);
sum = darrayCreate(3, darrayGetDimensions(dst));
if (!sum)
/* clear dst array with zeros, fill sum array with ones */
pixles = w*h*kd;
dst_ptr = darrayGetPr(dst);
- for (i = 0; i < pixles; i++)
+ for (i = 0; i < pixles; i++) {
+ darrayCheckPr(dst, dst_ptr, 1);
*dst_ptr++ = 0.0;
+ }
sum_ptr = darrayGetPr(sum);
- for (i = 0; i < pixles; i++)
+ for (i = 0; i < pixles; i++) {
+ darrayCheckPr(sum, sum_ptr, 1);
*sum_ptr++ = 1.0;
+ }
/* first copy all even planes */
for (z = 0; z < ks; z++) {
if ((z&1))
continue;
+ darrayCheckPr(src, darrayGetPr(src) + z*w*h, w*h);
+ darrayCheckPr(dst, darrayGetPr(dst) + (z/2)*w*h, w*h);
memcpy(darrayGetPr(dst) + (z/2)*w*h, darrayGetPr(src) + z*w*h, w*h*sizeof(double));
}
src_ptr = darrayGetPr(src) + z*w*h;
dst_ptr = darrayGetPr(dst) + (z/2)*w*h;
sum_ptr = darrayGetPr(sum) + (z/2)*w*h;
- flow_x_ptr = darrayGetPr(flow_x) + (z-1)*w*h;
- flow_y_ptr = darrayGetPr(flow_y) + (z-1)*w*h;
+ flow_x_ptr = darrayGetPr(flow) + (z-1)*w*h;
+ flow_y_ptr = flow_x_ptr + (ks-1)*w*h;
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
+ darrayCheckPr(flow, flow_x_ptr, 1);
xx = x + *flow_x_ptr++;
+ darrayCheckPr(flow, flow_y_ptr, 1);
yy = y + *flow_y_ptr++;
if (xx>=0 && yy>=0 && xx<w && yy<h) {
+ darrayCheckPr(dst, dst_ptr, 1);
+ darrayCheckPr(src, src_ptr + xx + yy*w, 1);
*dst_ptr += src_ptr[xx + yy*w];
+ darrayCheckPr(sum, sum_ptr, 1);
(*sum_ptr) += 1.0;
}
dst_ptr++;
pixles = w*h*kd;
dst_ptr = darrayGetPr(dst);
sum_ptr = darrayGetPr(sum);
- for (i = 0; i < pixles; i++)
+ for (i = 0; i < pixles; i++) {
+ darrayCheckPr(dst, dst_ptr, 1);
+ darrayCheckPr(sum, sum_ptr, 1);
*dst_ptr++ /= *sum_ptr++;
+ }
darrayDestroy(sum);
if (dst_move)
* A = 0.707
* B = 0.293
*/
-static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x, const darray_t *flow_y)
+static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow)
{
int x, y, z, xx, yy, xxx, yyy, w, h, ks, kd, k_;
double _x, _y;
ASSERT(w == darrayGetDimensions(dst)[0], "x-size of source and destination must be equal");
ASSERT(h == darrayGetDimensions(dst)[1], "y-size of source and destination must be equal");
ASSERT(ks == ((kd+1)>>1), "z-size of source must be half of the destination (+an extra plane if destination z-size is odd)");
- ASSERT_FLOW(w, h, kd, flow_x, flow_y, NULL, NULL);
+ ASSERT_FLOW(w, h, kd, flow, NULL);
sum = darrayCreate(3, darrayGetDimensions(dst));
if (!sum)
k_ = kd&(~1); /* even number of planes */
if (k_) {
/* copy first plane with weight 1.0 */
+ darrayCheckPr(dst, darrayGetPr(dst) + z*w*h, w*h);
+ darrayCheckPr(src, darrayGetPr(src) + (z/2)*w*h, w*h);
memcpy(darrayGetPr(dst) + z*w*h, darrayGetPr(src) + (z/2)*w*h, w*h*sizeof(double));
/* fill sum array with ones */
pixles = w*h;
sum2_ptr = sum1_ptr + w*h;
src1_ptr = darrayGetPr(src) + (z/2)*w*h;
src2_ptr = src1_ptr + w*h;
- flow_x1_ptr = darrayGetPr(flow_x) + z*w*h;
- flow_y1_ptr = darrayGetPr(flow_y) + z*w*h;
+ flow_x1_ptr = darrayGetPr(flow) + z*w*h;
+ flow_y1_ptr = flow_x1_ptr + (kd-1)*w*h;
flow_x2_ptr = flow_x1_ptr + w*h;
flow_y2_ptr = flow_y1_ptr + w*h;
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
/* even plane (weight A) */
+ darrayCheckPr(dst, dst2_ptr + x + y*w, 1);
+ darrayCheckPr(src, src2_ptr + x + y*w, 1);
dst2_ptr[x + y*w] += src2_ptr[x + y*w] * A;
+ darrayCheckPr(sum, sum2_ptr + x + y*w, 1);
sum2_ptr[x + y*w] += A;
/* get position via firt plane */
+ darrayCheckPr(flow, flow_x1_ptr, 1);
_x = *flow_x1_ptr++;
+ darrayCheckPr(flow, flow_y1_ptr, 1);
_y = *flow_y1_ptr++;
xx = x + round(_x);
yy = y + round(_y);
if (xx<0 || yy<0 || xx>=w || yy>=h)
continue;
/* odd plane (weight A)*/
+ darrayCheckPr(dst, dst1_ptr + xx + yy*w, 1);
+ darrayCheckPr(src, src1_ptr + x + y*w, 1);
dst1_ptr[xx + yy*w] += src1_ptr[x + y*w] * A;
+ darrayCheckPr(sum, sum1_ptr + xx + yy*w, 1);
sum1_ptr[xx + yy*w] += A;
/* follow position via second plane */
+ darrayCheckPr(flow, flow_x2_ptr + xx + yy*w, 1);
_x += flow_x2_ptr[xx + yy*w];
+ darrayCheckPr(flow, flow_y2_ptr + xx + yy*w, 1);
_y += flow_y2_ptr[xx + yy*w];
xxx = x + round(_x);
yyy = y + round(_y);
if (xxx<0 || yyy<0 || xxx>=w || yyy>=h)
continue;
/* even plane (weight B) */
+ darrayCheckPr(dst, dst2_ptr + xxx + yyy*w, 1);
+ darrayCheckPr(src, src1_ptr + x + y*w, 1);
dst2_ptr[xxx + yyy*w] += src1_ptr[x + y*w] * B;
+ darrayCheckPr(sum, sum2_ptr + xxx + yyy*w, 1);
sum2_ptr[xxx + yyy*w] += B;
/* odd plane (weigth B) */
+ darrayCheckPr(dst, dst1_ptr + xx + yy*w, 1);
+ darrayCheckPr(src, src2_ptr + xxx + yyy*w, 1);
dst1_ptr[xx + yy*w] += src2_ptr[xxx + yyy*w] * B;
+ darrayCheckPr(sum, sum1_ptr + xx + yy*w, 1);
sum1_ptr[xx + yy*w] += B;
}
}
dst1_ptr = darrayGetPr(dst) + z*w*h;
sum1_ptr = darrayGetPr(sum) + z*w*h;
src1_ptr = darrayGetPr(src) + (z/2)*w*h;
- flow_x1_ptr = darrayGetPr(flow_x) + z*w*h;
- flow_y1_ptr = darrayGetPr(flow_y) + z*w*h;
+ flow_x1_ptr = darrayGetPr(flow) + z*w*h;
+ flow_y1_ptr = flow_x1_ptr + (kd-1)*w*h;
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
/* get position via firt plane */
+ darrayCheckPr(flow, flow_x1_ptr, 1);
_x = *flow_x1_ptr++;
+ darrayCheckPr(flow, flow_y1_ptr, 1);
_y = *flow_y1_ptr++;
xx = x + round(_x);
yy = y + round(_y);
if (xx<0 || yy<0 || xx>=w || yy>=h)
continue;
/* border (odd) plane (weight 1.0)*/
+ darrayCheckPr(dst, dst1_ptr + xx + yy*w, 1);
+ darrayCheckPr(src, src1_ptr + x + y*w, 1);
dst1_ptr[xx + yy*w] += src1_ptr[x + y*w];
+ darrayCheckPr(sum, sum1_ptr + xx + yy*w, 1);
sum1_ptr[xx + yy*w] += 1.0;
}
}
/* copy extra plane, if destiantion size is odd */
if (z < kd) {
+ darrayCheckPr(dst, darrayGetPr(dst) + z*w*h, w*h);
+ darrayCheckPr(src, darrayGetPr(src) + (z/2)*w*h, w*h);
memcpy(darrayGetPr(dst) + z*w*h, darrayGetPr(src) + (z/2)*w*h, w*h*sizeof(double));
/* fill sum array with ones */
pixles = w*h;
sum1_ptr = darrayGetPr(sum);
for (i = 0; i < pixles; i++) {
/* only if there is a sum */
+ darrayCheckPr(sum, sum1_ptr, 1);
+ darrayCheckPr(dst, dst1_ptr, 1);
if (*sum1_ptr > 0.1)
*dst1_ptr++ /= *sum1_ptr;
else
*
* If the number of flow planes is odd, the extra plane is removed.
*/
-static int shrink_flow_z(darray_t *src_x, darray_t *src_y, darray_t *dst_x, darray_t *dst_y)
+int shrink_flow_z(darray_t *src, darray_t *dst)
{
int w, h, ks, kd, x, y, z, xx, yy;
double _x, _y;
double *src_x1_ptr, *src_x2_ptr, *src_y1_ptr, *src_y2_ptr, *dst_x_ptr, *dst_y_ptr;
- int dst_x_move = 0, dst_y_move = 0;
+ int dst_move = 0;
- ASSERT(darrayGetNumberOfDimensions(src_x) == 3, "src_x array must have 3 dimensions");
- ASSERT(darrayGetNumberOfDimensions(src_y) == 3, "src_y array must have 3 dimensions");
- w = darrayGetDimensions(src_x)[0];
- h = darrayGetDimensions(src_x)[1];
- ks = darrayGetDimensions(src_x)[2];
+ ASSERT(darrayGetNumberOfDimensions(src) == 4, "src array must have 4 dimensions");
+ w = darrayGetDimensions(src)[0];
+ h = darrayGetDimensions(src)[1];
+ ks = darrayGetDimensions(src)[2];
+ ASSERT(2 == darrayGetDimensions(src)[3], "number of cubes must be 2");
kd = ks>>1;
- ASSERT(w == darrayGetDimensions(src_y)[0], "x-size must be equal");
- ASSERT(h == darrayGetDimensions(src_y)[1], "y-size must be equal");
- ASSERT(ks == darrayGetDimensions(src_y)[2], "z-size must be equal");
-
- if (dst_x == NULL) {
- int dimensions[3] = {w, h, kd};
- dst_x = darrayCreate(3, dimensions);
- if (!dst_x)
- goto error;
- dst_x_move = 1;
- }
- if (dst_y == NULL) {
- int dimensions[3] = {w, h, kd};
+ if (dst == NULL) {
+ int dimensions[4] = {w, h, kd, 2};
- dst_y = darrayCreate(3, dimensions);
- if (!dst_y)
+ dst = darrayCreate(4, dimensions);
+ if (!dst)
goto error;
- dst_y_move = 1;
+ dst_move = 1;
}
- ASSERT(darrayGetNumberOfDimensions(dst_x) == 3, "flow_x array must have 3 dimensions");
- ASSERT(darrayGetNumberOfDimensions(dst_y) == 3, "flow_y array must have 3 dimensions");
- ASSERT(w == darrayGetDimensions(dst_x)[0], "x-size must be equal");
- ASSERT(h == darrayGetDimensions(dst_x)[1], "y-size must be equal");
- ASSERT(kd == darrayGetDimensions(dst_x)[2], "z-size of dest must be half of the planes + optional odd plane");
- ASSERT(w == darrayGetDimensions(dst_y)[0], "x-size must be equal");
- ASSERT(h == darrayGetDimensions(dst_y)[1], "y-size must be equal");
- ASSERT(kd == darrayGetDimensions(dst_y)[2], "z-size of dest must be half of the planes + optional odd plane");
+ ASSERT(darrayGetNumberOfDimensions(dst) == 4, "dst array must have 4 dimensions");
+ ASSERT(w == darrayGetDimensions(dst)[0], "x-size must be equal");
+ ASSERT(h == darrayGetDimensions(dst)[1], "y-size must be equal");
+ ASSERT(kd == darrayGetDimensions(dst)[2], "z-size of dest must be half of the planes + optional odd plane");
+ ASSERT(2 == darrayGetDimensions(dst)[3], "number of cubes must be 2");
for (z = 0; z < kd; z++) {
- src_x1_ptr = darrayGetPr(src_x) + (z*2)*w*h;
- src_y1_ptr = darrayGetPr(src_y) + (z*2)*w*h;
- dst_x_ptr = darrayGetPr(dst_x) + z*w*h;
- dst_y_ptr = darrayGetPr(dst_y) + z*w*h;
- src_x2_ptr = darrayGetPr(src_x) + (z*2+1)*w*h;
- src_y2_ptr = darrayGetPr(src_y) + (z*2+1)*w*h;
+ src_x1_ptr = darrayGetPr(src) + (z*2)*w*h;
+ src_y1_ptr = src_x1_ptr + w*h*ks;
+ src_x2_ptr = src_x1_ptr + w*h;
+ src_y2_ptr = src_y1_ptr + w*h;
+ dst_x_ptr = darrayGetPr(dst) + z*w*h;
+ dst_y_ptr = dst_x_ptr + w*h*kd;
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
/* get position via first plane */
+ darrayCheckPr(src, src_x1_ptr, 1);
_x = *src_x1_ptr++;
+ darrayCheckPr(src, src_y1_ptr, 1);
_y = *src_y1_ptr++;
xx = x + round(_x);
yy = y + round(_y);
if (xx<0 || yy<0 || xx>=w || yy>=h) {
/* use an outside value when outside in the first plane */
+ darrayCheckPr(dst, dst_x_ptr, 1);
*dst_x_ptr++ = _x;
+ darrayCheckPr(dst, dst_y_ptr, 1);
*dst_y_ptr++ = _y;
continue;
}
/* follow position via second plane */
+ darrayCheckPr(src, src_x2_ptr + xx + yy*w, 1);
_x += src_x2_ptr[xx + yy*w];
+ darrayCheckPr(src, src_y2_ptr + xx + yy*w, 1);
_y += src_y2_ptr[xx + yy*w];
/* store displacement result */
+ darrayCheckPr(dst, dst_x_ptr, 1);
*dst_x_ptr++ = _x;
+ darrayCheckPr(dst, dst_y_ptr, 1);
*dst_y_ptr++ = _y;
}
}
}
- if (dst_x_move)
- darrayMove(dst_x, src_x);
- if (dst_y_move)
- darrayMove(dst_y, src_y);
+ if (dst_move)
+ darrayMove(dst, src);
return 0;
error:
- if (dst_x_move)
- darrayDestroy(dst_x);
- if (dst_y_move)
- darrayDestroy(dst_y);
+ if (dst_move)
+ darrayDestroy(dst);
return -1;
}
* If a dimension of a source matrix is odd, the end of the dimension is just
* copied and not scaled.
*
- * The points in a plane are shrinked from 2x2x2 to 1x1x1 by calculating the mean
+ * The points in a plane are shrinked from 2x2 to 1x1 by calculating the mean
* of all eight values.
+ *
+ * If the array has more than two dimensions, all planes in all dimensions
+ * are processed.
*/
-static int shrink_xy(darray_t *src, darray_t *dst)
+int shrink_xy(darray_t *src, darray_t *dst)
{
- int x, xx, y, z, ws, wd, w_, hs, hd, h_, k;
+ int x, xx, y, z, ws, wd, w_, hs, hd, h_, k, dims, i;
double *s1, *s2, *d;
int dst_move = 0;
-
- ASSERT(darrayGetNumberOfDimensions(src) == 3, "Source array must have 3 dimensions");
+
+ dims = darrayGetNumberOfDimensions(src);
+ ASSERT(dims >= 2, "Source array must have at least 2 dimensions");
ws = darrayGetDimensions(src)[0];
w_ = ws>>1;
hs = darrayGetDimensions(src)[1];
h_ = hs>>1;
- k = darrayGetDimensions(src)[2];
+ k = 1;
+ for (i = 2; i < dims; i++)
+ k *= darrayGetDimensions(src)[i];
/* create dst, if not given */
if (!dst) {
- int dimensions[3] = {(ws+1)>>1, (hs+1)>>1, k};
-
- dst = darrayCreate(3, dimensions);
+ int dimensions[dims];
+
+ dimensions[0] = (ws+1)>>1;
+ dimensions[1] = (hs+1)>>1;
+ for (i = 2; i < dims; i++)
+ dimensions[i] = darrayGetDimensions(src)[i];
+
+ dst = darrayCreate(dims, dimensions);
if (!dst)
goto error;
dst_move = 1;
}
- ASSERT(darrayGetNumberOfDimensions(dst) == 3, "Destination array must have 3 dimensions");
+ ASSERT(darrayGetNumberOfDimensions(dst) == dims, "Destination array must have equal dimensions than source array");
wd = darrayGetDimensions(dst)[0];
hd = darrayGetDimensions(dst)[1];
ASSERT((ws&1) || ws == wd*2, "x-size of source must be twice the size of destination (this is even source size)");
ASSERT(!(ws&1) || ws-1 == (wd-1)*2, "x-size of source-1 must be twice the size of destination-1 (this odd source size)");
ASSERT((hs&1) || hs == hd*2, "y-size of source must be twice the size of destination (this is even source size)");
ASSERT(!(hs&1) || hs-1 == (hd-1)*2, "y-size of source-1 must be twice the size of destination-1 (this odd source size)");
- ASSERT(k == darrayGetDimensions(dst)[2], "z-size must be equal");
+ for (i = 2; i < dims; i++)
+ ASSERT(darrayGetDimensions(src)[i] == darrayGetDimensions(dst)[i], "dimension 3 and greater must be equal");
for (z = 0; z < k; z++) {
for (y = 0; y < h_; y++) {
s1 = darrayGetPr(src) + y*2*ws + z*ws*hs;
s2 = s1 + ws;
d = darrayGetPr(dst) + y*wd + z*wd*hd;
- for (x = 0, xx = 0; x < w_; x++, xx+=2)
+ for (x = 0, xx = 0; x < w_; x++, xx+=2) {
+ darrayCheckPr(dst, d + x, 1);
+ darrayCheckPr(src, s1 + xx, 2);
+ darrayCheckPr(src, s2 + xx, 2);
d[x] = (s1[xx] + s1[xx+1] + s2[xx] + s2[xx+1]) * 0.25;
+ }
/* mean of 1x2x1 */
- if (x < wd)
+ if (x < wd) {
+ darrayCheckPr(dst, d + x, 1);
+ darrayCheckPr(src, s1 + xx, 1);
+ darrayCheckPr(src, s2 + xx, 1);
d[x] = (s1[xx] + s2[xx]) * 0.5;
+ }
}
if (y < hd) {
/* mean of 2x1x1 */
s1 = darrayGetPr(src) + y*2*ws + z*ws*hs;
d = darrayGetPr(dst) + y*wd + z*wd*hd;
- for (x = 0, xx = 0; x < w_; x++, xx+=2)
+ for (x = 0, xx = 0; x < w_; x++, xx+=2) {
+ darrayCheckPr(dst, d + x, 1);
+ darrayCheckPr(src, s1 + xx, 2);
d[x] = (s1[xx] + s1[xx+1]) * 0.5;
+ }
/* mean of 1x1x1 */
- if (x < wd)
+ if (x < wd) {
+ darrayCheckPr(dst, d + x, 1);
+ darrayCheckPr(src, s1 + xx, 1);
d[x] = s1[xx];
+ }
}
}
d = darrayGetPr(dst) + (y+h_)*wd + z*wd*hd;
*d++ = *s; /* x=0 plane */
for (x = 0; x < w_-1; x++) {
+ darrayCheckPr(dst, d, 1);
+ darrayCheckPr(src, s, 2);
*d++ = A*s[0] + B*s[1];
*d++ = A*s[1] + B*s[0];
s++;
}
+ darrayCheckPr(dst, d, 1);
+ darrayCheckPr(src, s, 1);
*d++ = *s++; /* x=w-1 plane */
/* extra plane (odd fine grid) */
- if ((wd&1))
+ if ((wd&1)) {
+ darrayCheckPr(dst, d, 1);
+ darrayCheckPr(src, s, 1);
*d = *s;
+ }
}
}
/* y=h-1 plane */
s = darrayGetPr(dst) + h_*wd + z*wd*hd;
d = darrayGetPr(dst) + z*wd*hd;
- for (x = 0; x < wd; x++)
+ for (x = 0; x < wd; x++) {
+ darrayCheckPr(dst, d, 1);
+ darrayCheckPr(dst, s, 1);
*d++ = *s++; /* y=0 plane */
+ }
for (y = 0; y < h_-1; y++) {
s1 = darrayGetPr(dst) + (y+h_)*wd + z*wd*hd;
s2 = s1 + wd;
d1 = darrayGetPr(dst) + (y*2+1)*wd + z*wd*hd;
d2 = d1 + wd;
for (x = 0; x < wd; x++) {
+ darrayCheckPr(dst, d1, 1);
+ darrayCheckPr(dst, d2, 1);
+ darrayCheckPr(dst, s1, 1);
+ darrayCheckPr(dst, s2, 1);
*d1++ = A * *s1 + B * *s2;
*d2++ = A * *s2 + B * *s1;
s1++;
*
* If an element is greater theshold, the element is set to gt, otherwise le.
*/
-static void quantize(const darray_t *array, double threshold, double gt, double le)
+void quantize(const darray_t *array, double threshold, double gt, double le)
{
- int n = 1, i;
+ int n = 1, dims, i;
double *ptr;
- ASSERT(darrayGetNumberOfDimensions(array) == 3, "array must have 3 dimensions");
- n *= darrayGetDimensions(array)[0];
- n *= darrayGetDimensions(array)[1];
- n *= darrayGetDimensions(array)[2];
+ dims = darrayGetNumberOfDimensions(array);
+ for (i = 0; i < dims; i++)
+ n *= darrayGetDimensions(array)[i];
ptr = darrayGetPr(array);
for (i = 0; i < n; i++, ptr++) {
+ darrayCheckPr(array, ptr, 1);
if (*ptr > threshold)
*ptr = gt;
else
/* multiply the array content with a value
*/
-static void multiply(const darray_t *array, double factor)
+void multiply(const darray_t *array, double factor)
{
- int n = 1, i;
+ int n = 1, dims, i;
double *ptr;
- ASSERT(darrayGetNumberOfDimensions(array) == 3, "array must have 3 dimensions");
- n *= darrayGetDimensions(array)[0];
- n *= darrayGetDimensions(array)[1];
- n *= darrayGetDimensions(array)[2];
+ dims = darrayGetNumberOfDimensions(array);
+ for (i = 0; i < dims; i++)
+ n *= darrayGetDimensions(array)[i];
ptr = darrayGetPr(array);
- for (i = 0; i < n; i++, ptr++)
+ for (i = 0; i < n; i++, ptr++) {
+ darrayCheckPr(array, ptr, 1);
*ptr *= factor;
+ }
}
/* create an array of neighbors and use flow vectors
+ *
+ * the elements of the array consist of:
+ * int num
+ * int index[neighbors]
+ * double weight[neighbors]
+ * num specifies the number of neighbors used.
+ * neighbors is fixed for all entries in the array and depends on 2d or 3d
+ * image. since the number of neighbors depends on the dimensions of an image,
+ * they are not fixed. a structure cannot be used in this case.
+ *
* flow_x and flow_x are the vectors that point to neighbors of the next z plane
* flow_x and flow_x are the vectors that point to neighbors of the previous z plane
*/
-static unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy)
+unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow, const darray_t *flow_i)
{
int x, y, z, xx, yy, zz, i, j, n;
double _x, _y;
int index;
double *flow_x_ptr, *flow_y_ptr, *flow_ix_ptr, *flow_iy_ptr;
unsigned char *neighbors, *ret_neighbors;
+ int nb_num, index_offset;
unsigned long nb_size;
- int nb_num;
int *i_index;
- ASSERT_FLOW(w, h, k, flow_x, flow_y, flow_ix, flow_iy);
- flow_x_ptr = darrayGetPr(flow_x);
- flow_y_ptr = darrayGetPr(flow_y);
- flow_ix_ptr = darrayGetPr(flow_ix);
- flow_iy_ptr = darrayGetPr(flow_iy);
+ ASSERT_FLOW(w, h, k, flow, flow_i);
+ flow_x_ptr = darrayGetPr(flow);
+ flow_y_ptr = flow_x_ptr + w*h*(k-1);
+ flow_ix_ptr = darrayGetPr(flow_i);
+ flow_iy_ptr = flow_ix_ptr + w*h*(k-1);
nb_num = (k==1) ? NEIGHBORS_2D : NEIGHBORS_3D;
- nb_size = (k==1) ? sizeof(struct neighbors_2d) : sizeof(struct neighbors_3d);
+ nb_size = sizeof(int) + (sizeof(int) + sizeof(double)) * nb_num;
+ index_offset = sizeof(int);
ret_neighbors = neighbors = calloc(w*h*k, nb_size);
// printf("alloc neighbor grid=%lu Mbytes\n", nb_size*w*h*k/1024/1024);
if (!ret_neighbors) {
- printf("failed to allocate neighbor grid\n");
+ printf("failed to allocate neighbor grid (%lu megabytes)\n", nb_size * w*h*k / 1024 / 1024 );
return NULL;
}
for (z = 0; z < k; z++) {
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
- if (k==1)
- i_index = ((struct neighbors_2d *)neighbors)->index;
- else
- i_index = ((struct neighbors_3d *)neighbors)->index;
+ i_index = (int *)(neighbors+index_offset);
for (i = 0, n = 0; i < nb_num; i++) {
/* xx, yy, zz is the neighbor of color source
* The first neighbor is always the point itself!
if (xx >= w || yy >= h || zz >= k)
continue;
if (zz > z) {
+ darrayCheckPr(flow, flow_x_ptr + xx + yy*w + (zz-1)*w*h, 1);
_x = flow_x_ptr[xx + yy*w + (zz-1)*w*h];
+ darrayCheckPr(flow, flow_y_ptr + xx + yy*w + (zz-1)*w*h, 1);
_y = flow_y_ptr[xx + yy*w + (zz-1)*w*h];
xx = xx + round(_x);
yy = yy + round(_y);
} else if (zz < z) {
+ darrayCheckPr(flow_i, flow_ix_ptr + xx + yy*w + zz*w*h, 1);
_x = flow_ix_ptr[xx + yy*w + zz*w*h];
+ darrayCheckPr(flow_i, flow_iy_ptr + xx + yy*w + zz*w*h, 1);
_y = flow_iy_ptr[xx + yy*w + zz*w*h];
xx = xx + round(_x);
yy = yy + round(_y);
continue;
i_index[n++] = index;
}
- if (k == 1)
- ((struct neighbors_2d *)neighbors)->num = n;
- else
- ((struct neighbors_3d *)neighbors)->num = n;
+ *((int *)neighbors) = n;
neighbors += nb_size;
}
}
}
/* weighten neighbors */
-static void weighten(const darray_t *luminance, unsigned char *neighbors)
+void weighten(const darray_t *luminance, unsigned char *neighbors)
{
int x, y, z, i, w, h, k, num;
double val, nval, mean, variance, sigma, sum, weight, min;
double sqdiff[NEIGHBORS_3D];
double log0_01 = log(0.01);
double *lum_ptr;
+ int nb_num, index_offset, weight_offset;
unsigned long nb_size;
int *i_index; double *w_index;
k = darrayGetDimensions(luminance)[2];
lum_ptr = darrayGetPr(luminance);
- nb_size = (k==1) ? sizeof(struct neighbors_2d) : sizeof(struct neighbors_3d);
+ nb_num = (k==1) ? NEIGHBORS_2D : NEIGHBORS_3D;
+ nb_size = sizeof(int) + (sizeof(int) + sizeof(double)) * nb_num;
+ index_offset = sizeof(int);
+ weight_offset = sizeof(int) + sizeof(int) * nb_num;
for (z = 0; z < k; z++) {
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++, neighbors += nb_size) {
- if (k==1) {
- num = ((struct neighbors_2d *)neighbors)->num;
- i_index = ((struct neighbors_2d *)neighbors)->index;
- w_index = ((struct neighbors_2d *)neighbors)->weight;
- } else {
- num = ((struct neighbors_3d *)neighbors)->num;
- i_index = ((struct neighbors_3d *)neighbors)->index;
- w_index = ((struct neighbors_3d *)neighbors)->weight;
- }
+ num = *((int *)neighbors);
+ i_index = (int *)(neighbors+index_offset);
+ w_index = (double *)(neighbors+weight_offset);
/* Be sure that we have at least 3 neighbors.
* This applies for every corner in a plane.
* The first point in the neighbor list iss the point itself.
variance = 0;
for (i = 0; i < num; i++) {
/* value of neighbor pixles, including current pixle */
+ darrayCheckPr(luminance, lum_ptr + i_index[i], 1);
nval = lum_ptr[i_index[i]];
/* squared diffs between current pixle and each neighbor */
sqdiff[i] = (nval-val) * (nval-val);
int x, y, z, w, h, k, i, num;
double *value_ptr, *value_line_ptr, *mark_line_ptr, *init_line_ptr;
double val;
+ int nb_num, index_offset, weight_offset;
unsigned long nb_size;
int *i_index; double *w_index;
ASSERT(k == darrayGetDimensions(mark)[2], "z-sizes of mark array is different");
value_ptr = darrayGetPr(value);
- nb_size = (k==1) ? sizeof(struct neighbors_2d) : sizeof(struct neighbors_3d);
+ nb_num = (k==1) ? NEIGHBORS_2D : NEIGHBORS_3D;
+ nb_size = sizeof(int) + (sizeof(int) + sizeof(double)) * nb_num;
+ index_offset = sizeof(int);
+ weight_offset = sizeof(int) + sizeof(int) * nb_num;
for (z = 0; z < k; z++) {
for (y = 0; y < h; y++) {
mark_line_ptr = darrayGetPr(mark) + y*w + z*w*h;
init_line_ptr = darrayGetPr(init) + y*w + z*w*h;
for (x = 0; x < w; x++, neighbors += nb_size) {
- if (k==1) {
- num = ((struct neighbors_2d *)neighbors)->num;
- i_index = ((struct neighbors_2d *)neighbors)->index;
- w_index = ((struct neighbors_2d *)neighbors)->weight;
- } else {
- num = ((struct neighbors_3d *)neighbors)->num;
- i_index = ((struct neighbors_3d *)neighbors)->index;
- w_index = ((struct neighbors_3d *)neighbors)->weight;
- }
+ num = *((int *)neighbors);
+ i_index = (int *)(neighbors+index_offset);
+ w_index = (double *)(neighbors+weight_offset);
/* if pixle is marked */
+ darrayCheckPr(mark, mark_line_ptr + x, 1);
if (mark_line_ptr[x] < 0.5) {
/* always use initial color on marked plane, otherwise do not smooth */
- if (level == 0)
+ if (level == 0) {
+ darrayCheckPr(value, value_line_ptr + x, 1);
+ darrayCheckPr(init, init_line_ptr + x, 1);
value_line_ptr[x] = init_line_ptr[x];
+ }
continue;
}
/* Calculate value */
val = 0;
- for (i = 1; i < num; i++) /* sum value without point itself */
+ for (i = 1; i < num; i++) { /* sum value without point itself */
+ darrayCheckPr(value, value_ptr + i_index[i], 1);
val += w_index[i] * value_ptr[i_index[i]];
- if (level == 0)
+ }
+ if (level == 0) {
+ darrayCheckPr(value, value_line_ptr + x, 1);
+ darrayCheckPr(init, init_line_ptr + x, 1);
value_line_ptr[x] = (init_line_ptr[x] - val) / w_index[0];
- else
+ } else {
+ darrayCheckPr(value, value_line_ptr + x, 1);
value_line_ptr[x] = (-val) / w_index[0];
+ }
}
}
}
}
+#ifdef RESIDUAL
/* Calculate residual
*/
static double residual(int level, darray_t *value, const darray_t *mark, const darray_t *init, unsigned char *neighbors)
int x, y, z, w, h, k, i, num;
double *value_ptr, *mark_line_ptr, *init_line_ptr;
double val, residual = 0;
+ int nb_num, index_offset, weight_offset;
unsigned long nb_size;
int *i_index; double *w_index;
ASSERT(k == darrayGetDimensions(mark)[2], "z-sizes of mark array is different");
value_ptr = darrayGetPr(value);
- nb_size = (k==1) ? sizeof(struct neighbors_2d) : sizeof(struct neighbors_3d);
+ nb_num = (k==1) ? NEIGHBORS_2D : NEIGHBORS_3D;
+ nb_size = sizeof(int) + (sizeof(int) + sizeof(double)) * nb_num;
+ index_offset = sizeof(int);
+ weight_offset = sizeof(int) + sizeof(int) * nb_num;
for (z = 0; z < k; z++) {
for (y = 0; y < h; y++) {
init_line_ptr = darrayGetPr(init) + y*w + z*w*h;
mark_line_ptr = darrayGetPr(mark) + y*w + z*w*h;
for (x = 0; x < w; x++, neighbors += nb_size) {
- if (k==1) {
- num = ((struct neighbors_2d *)neighbors)->num;
- i_index = ((struct neighbors_2d *)neighbors)->index;
- w_index = ((struct neighbors_2d *)neighbors)->weight;
- } else {
- num = ((struct neighbors_3d *)neighbors)->num;
- i_index = ((struct neighbors_3d *)neighbors)->index;
- w_index = ((struct neighbors_3d *)neighbors)->weight;
- }
+ num = *((int *)neighbors);
+ i_index = (int *)(neighbors+index_offset);
+ w_index = (double *)(neighbors+weight_offset);
/* if pixle is marked, no residual */
+ darrayCheckPr(mark, mark_line_ptr + x, 1);
if (mark_line_ptr[x] < 0.5)
continue;
/* Calculate value */
val = 0;
- for (i = 0; i < num; i++) /* sum value including point itself */
+ for (i = 0; i < num; i++) { /* sum value including point itself */
+ darrayCheckPr(value, value_ptr + i_index[i], 1);
val += w_index[i] * value_ptr[i_index[i]];
- if (level == 0)
+ }
+ if (level == 0) {
+ darrayCheckPr(init, init_line_ptr + x, 1);
val = init_line_ptr[x] - val;
- else
+ } else
val = -val;
residual += val * val;
}
return residual;
}
+#endif
/* Restrict
*
* different depth (z dimension), a temporary array (temp) is used to
* perform the intermediate step of scaling x,y,z.
*
- * Note that the flow vectors must be scaled xy already!
+ * Note that the flow vectors must be scaled xy already (if x and y is shrinked
+ * before z)!
*/
-static int restrict(const darray_t *value, darray_t *nvalue, const darray_t *flow_x, const darray_t *flow_y)
+static int restrict(const darray_t *value, darray_t *nvalue, const darray_t *flow)
{
int ks, wd, hd, kd;
darray_t *temp = NULL;
kd = darrayGetDimensions(nvalue)[2];
ks = darrayGetDimensions(value)[2];
if (ks != kd)
- ASSERT(flow_x && flow_y, "when scaling z, flow vectors must be given");
+ ASSERT(flow, "when scaling z, flow vectors must be given");
if (ks != kd) {
int dimensions[3] = {wd, hd, ks};
rc = shrink_xy((darray_t *)value, temp);
if (rc < 0)
goto error;
- rc = shrink_z(temp, nvalue, flow_x, flow_y);
+ rc = shrink_z(temp, nvalue, flow);
if (rc < 0)
goto error;
} else {
* different depth (z dimension), a temporary array (temp3) is used to
* perform the intermediate step of scaling x,y,z.
*
- * Note that the flow vectors must be scaled xy already!
+ * Note that the flow vectors must be scaled xy already (if x and y is shrinked
+ * before z)!
*
* First the value is shrinked to the size of the value of the next level.
* The error is calculated by substracting the values of the next level
* substracted from the value. So the changes in the next level's value are
* prolonged to the value.
*/
-static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow_x, const darray_t *flow_y)
+static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow)
{
int ws, hs, ks, wd, hd, kd;
int i, pixles;
hd = darrayGetDimensions(value)[1];
kd = darrayGetDimensions(value)[2];
if (ks != kd)
- ASSERT(flow_x && flow_y, "when scaling z, flow vectors must be given");
+ ASSERT(flow, "when scaling z, flow vectors must be given");
/* generate temporary arrays */
temp1 = darrayCreate(3, darrayGetDimensions(value));
rc = shrink_xy(value, temp3);
if (rc < 0)
goto error;
- rc = shrink_z(temp3, temp2, flow_x, flow_y);
+ rc = shrink_z(temp3, temp2, flow);
if (rc < 0)
goto error;
} else {
temp2_ptr = darrayGetPr(temp2);
nvalue_ptr = darrayGetPr(nvalue);
pixles = ws*hs*ks;
- for (i = 0; i < pixles; i++)
+ for (i = 0; i < pixles; i++) {
+ darrayCheckPr(temp2, temp2_ptr, 1);
+ darrayCheckPr(nvalue, nvalue_ptr, 1);
*temp2_ptr++ -= *nvalue_ptr++;
+ }
/* Expand (temp2 -> temp1) */
if (ks != kd) {
- rc = expand_z(temp2, temp3, flow_x, flow_y);
+ rc = expand_z(temp2, temp3, flow);
if (rc < 0)
goto error;
expand_xy(temp3, temp1);
temp1_ptr = darrayGetPr(temp1);
value_ptr = darrayGetPr(value);
pixles = wd*hd*kd;
- for (i = 0; i < pixles; i++)
+ for (i = 0; i < pixles; i++) {
+ darrayCheckPr(temp1, temp1_ptr, 1);
+ darrayCheckPr(value, value_ptr, 1);
*value_ptr++ -= *temp1_ptr++;
+ }
rc = 0;
for (z = 0; z < ks; z++) {
for (y = 0; y < hs; y++) {
for (x = 0; x < ws; x++) {
+ darrayCheckPr(dst, dst_ptr + (x+xoff) + (y+yoff)*wd + (z+zoff)*wd*hd, 1);
+ darrayCheckPr(src, src_ptr + x + y*ws + z*ws*hs, 1);
dst_ptr[(x+xoff) + (y+yoff)*wd + (z+zoff)*wd*hd] = src_ptr[x + y*ws + z*ws*hs];
}
}
* The flow vectors are drawn to the upper left corner of the dst, which
* may be larger in size than the flow arrays.
*/
-void draw_flow(darray_t *dst, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy)
+void draw_flow(darray_t *dst, const darray_t *flow, const darray_t *flow_i)
{
int wd, hd, ws, hs, ks, x, y, z;
double *flow_x_ptr, *flow_y_ptr, *flow_ix_ptr, *flow_iy_ptr;
wd = darrayGetDimensions(dst)[0];
hd = darrayGetDimensions(dst)[1];
- ws = darrayGetDimensions(flow_x)[0];
- hs = darrayGetDimensions(flow_x)[1];
- ks = darrayGetDimensions(flow_x)[2];
- flow_x_ptr = darrayGetPr(flow_x);
- flow_y_ptr = darrayGetPr(flow_y);
- flow_ix_ptr = darrayGetPr(flow_ix);
- flow_iy_ptr = darrayGetPr(flow_iy);
+ ws = darrayGetDimensions(flow)[0];
+ hs = darrayGetDimensions(flow)[1];
+ ks = darrayGetDimensions(flow)[2];
+ flow_x_ptr = darrayGetPr(flow);
+ flow_y_ptr = flow_x_ptr + ws*hs*ks;
+ flow_ix_ptr = darrayGetPr(flow_i);
+ flow_iy_ptr = flow_ix_ptr + ws*hs*ks;
for (z = 0; z <= ks; z++) {
for (y = step/2; y < hs-step/2+1; y+=step) {
/* Solve multigrid
*
- * Note that the flow vectors must be scaled xy already!
+ * Note that the flow vectors must be scaled xy already (if x and y is shrinked
+ * before z)!
+ *
*/
-static int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows_xy_x, darray_t **flows_xy_y, int scalexyz)
+int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows, int scalexyz, double target_residual_change)
{
+#ifdef RESIDUAL
double curr_residual, prev_residual = 10000000000000.0;
+#endif
int cycle, level, i;
int rc = -1;
for (cycle = 0; cycle < ncycle; cycle++) {
/* go down */
+#ifdef RESIDUAL
curr_residual = 0;
+#endif
for (level = 0; level < nlevel-1; level++) {
for (i = 0; i < iters; i++)
smooth(level, values[level], marks[level], init, nb_list[level]);
- curr_residual += residual(level, values[level], marks[level], init, nb_list[level]);
+#ifdef RESIDUAL
+ if (target_residual_change)
+ curr_residual += residual(level, values[level], marks[level], init, nb_list[level]);
+#endif
if (level < scalexyz)
- rc = restrict(values[level], values[level+1], flows_xy_x[level], flows_xy_y[level]);
+ rc = restrict(values[level], values[level+1], flows[level]);
else
- rc = restrict(values[level], values[level+1], NULL, NULL);
+ rc = restrict(values[level], values[level+1], NULL);
if (rc < 0)
goto error;
}
/* go up */
for (level = nlevel-2; level >= 0; level--) {
if (level < scalexyz)
- rc = prolong(values[level], values[level+1], flows_xy_x[level], flows_xy_y[level]);
+ rc = prolong(values[level], values[level+1], flows[level]);
else
- rc = prolong(values[level], values[level+1], NULL, NULL);
+ rc = prolong(values[level], values[level+1], NULL);
if (rc < 0)
goto error;
for (i = 0; i < iters; i++)
}
/* end if residual raises again */
-// if (prev_residual < curr_residual)
-// break;
- curr_residual = prev_residual;
+#ifdef RESIDUAL
+//printf("residual %.8f (change=%%%.4f)\n", curr_residual, (prev_residual-curr_residual)/prev_residual*100);
+ if (target_residual_change)
+ printf(" %.1f%%", (prev_residual-curr_residual)/prev_residual*100); fflush(stdout);
+ if (curr_residual == 0.0) {
+ printf(" 0"); fflush(stdout);
+ break;
+ }
+ if (target_residual_change && (prev_residual-curr_residual)/prev_residual < target_residual_change)
+ break;
+ prev_residual = curr_residual;
+#endif
}
return 0;
return rc;
}
-/* colorize image array using image_mark
- *
- * image:
- * first dimansion = width
- * second dimension = height
- * thrid dimension = components (first component is luminance, other components are features, like u, v)
- * fourth dimension = time
- *
- * image_mark:
- * first dimansion = width
- * second dimension = height
- * fourth dimension = time
- */
-darray_t *colorize(const darray_t *image, const darray_t *image_mark, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy, int inner_iter, int outer_iter, int scalexyz)
-{
- unsigned char **nb_list = NULL;
- darray_t **values = NULL, **marks = NULL, *init = NULL, **flows_x = NULL, **flows_y = NULL, **flows_ix = NULL, **flows_iy = NULL, **flows_xy_x = NULL, **flows_xy_y = NULL;
- darray_t *result = NULL, *luminance = NULL;
- int dimensions[4];
- int w, h, k, ww, hh, kk, last_kk, comp, d, max_depth, temp, z, c;
- int rc = -1;
-
- ASSERT(darrayGetNumberOfDimensions(image) == 4, "image array must have 4 dimensions");
- w = darrayGetDimensions(image)[0];
- h = darrayGetDimensions(image)[1];
- comp = darrayGetDimensions(image)[2];
- k = darrayGetDimensions(image)[3];
- ASSERT(comp > 0, "image array requires at least one component (Y-component)");
-
- ASSERT(darrayGetNumberOfDimensions(image_mark) == 3, "image_mark array must have 3 dimensions");
- ASSERT(w == darrayGetDimensions(image_mark)[0], "image_mark and image must both have same x-size");
- ASSERT(h == darrayGetDimensions(image_mark)[1], "image_mark and image must both have same y-size");
- ASSERT(k == darrayGetDimensions(image_mark)[2], "image_mark and image must both have same z-size");
-
- /* Calculate maximum depth
- * Scale down dimensions until at least 2x2 remains
- */
- temp = min(w, h);
- for (max_depth = 0; temp > 2; max_depth++)
- temp=(temp+1)>>1;
- ASSERT(max_depth >= 2, "image too small");
-
- nb_list = calloc(max_depth, sizeof(struct neighbors *));
- if (!nb_list)
- goto error;
- marks = calloc(max_depth, sizeof(darray_t *));
- if (!marks)
- goto error;
- ASSERT(marks, "no memory");
- values = calloc(max_depth, sizeof(darray_t *));
- if (!values)
- goto error;
- ASSERT(values, "no memory");
- flows_x = calloc(max_depth, sizeof(darray_t *));
- if (!flows_x)
- goto error;
- ASSERT(flows_x, "no memory");
- flows_y = calloc(max_depth, sizeof(darray_t *));
- if (!flows_y)
- goto error;
- ASSERT(flows_y, "no memory");
- flows_ix = calloc(max_depth, sizeof(darray_t *));
- if (!flows_ix)
- goto error;
- ASSERT(flows_ix, "no memory");
- flows_iy = calloc(max_depth, sizeof(darray_t *));
- if (!flows_iy)
- goto error;
- ASSERT(flows_iy, "no memory");
- if (scalexyz) {
- /* the flows_xy are scaled in xy in the current level */
- flows_xy_x = calloc(max_depth, sizeof(darray_t *));
- if (!flows_xy_x)
- goto error;
- flows_xy_y = calloc(max_depth, sizeof(darray_t *));
- if (!flows_xy_y)
- goto error;
- }
- dimensions[0] = w;
- dimensions[1] = h;
- dimensions[2] = k;
- luminance = darrayCreate(3, dimensions);
- if (!luminance)
- goto error;
- /* generate sparse multigrid */
- last_kk = kk = k;
- printf("-> prepare sparse grid");
- for (d = 0, ww=w, hh=h; d < max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
- /* create temporary arrays
- * copy fist array to level 0 or shrink previous level
- */
- printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
- dimensions[0] = ww;
- dimensions[1] = hh;
- dimensions[2] = last_kk;
- marks[d] = darrayCreate(3, dimensions);
- if (!marks[d])
- goto error;
- dimensions[0] = ww;
- dimensions[1] = hh;
- dimensions[2] = last_kk-1;
- flows_x[d] = darrayCreate(3, dimensions);
- if (!flows_x[d])
- goto error;
- flows_y[d] = darrayCreate(3, dimensions);
- if (!flows_y[d])
- goto error;
- flows_ix[d] = darrayCreate(3, dimensions);
- if (!flows_ix[d])
- goto error;
- flows_iy[d] = darrayCreate(3, dimensions);
- if (!flows_iy[d])
- goto error;
- if (d < scalexyz) {
- dimensions[0] = (ww+1)>>1;
- dimensions[1] = (hh+1)>>1;
- dimensions[2] = kk-1;
- flows_xy_x[d] = darrayCreate(3, dimensions);
- if (!flows_xy_x[d])
- goto error;
- flows_xy_y[d] = darrayCreate(3, dimensions);
- if (!flows_xy_y[d])
- goto error;
- }
- if (d == 0) {
- for (z = 0; z < kk; z++)
- memcpy(darrayGetPr(luminance) + z*w*h, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
- if (flow_x && flow_y && flow_ix && flow_iy) {
- memcpy(darrayGetPr(flows_x[d]), darrayGetPr(flow_x), w*h*(kk-1)*sizeof(double));
- memcpy(darrayGetPr(flows_y[d]), darrayGetPr(flow_y), w*h*(kk-1)*sizeof(double));
- memcpy(darrayGetPr(flows_ix[d]), darrayGetPr(flow_ix), w*h*(kk-1)*sizeof(double));
- memcpy(darrayGetPr(flows_iy[d]), darrayGetPr(flow_iy), w*h*(kk-1)*sizeof(double));
- }
- memcpy(darrayGetPr(marks[d]), darrayGetPr(image_mark), w*h*kk*sizeof(double));
- /* invert, so 0 = marked */
- quantize(marks[d], 0.5, 0, 1);
- } else {
- rc = shrink_xy(luminance, NULL);
- if (rc < 0)
- goto error;
- rc = shrink_xy(marks[d-1], marks[d]);
- if (rc < 0)
- goto error;
- rc = shrink_xy(flows_x[d-1], flows_x[d]);
- if (rc < 0)
- goto error;
- rc = shrink_xy(flows_y[d-1], flows_y[d]);
- if (rc < 0)
- goto error;
- rc = shrink_xy(flows_ix[d-1], flows_ix[d]);
- if (rc < 0)
- goto error;
- rc = shrink_xy(flows_iy[d-1], flows_iy[d]);
- if (rc < 0)
- goto error;
- /* shrink length of flow vector on smaller grid.
- * This is not performed when comparing with
- * the original algorithm. (it is a bug there)
- */
- multiply(flows_x[d], 0.5);
- multiply(flows_y[d], 0.5);
- multiply(flows_ix[d], 0.5);
- multiply(flows_iy[d], 0.5);
- if (d-1 < scalexyz) { /* process if previous d shall be shrinked in z direction */
- /* shrink time */
- rc = shrink_z(luminance, NULL, flows_x[d], flows_y[d]);
- if (rc < 0)
- goto error;
- rc = shrink_z(marks[d], NULL, flows_x[d], flows_y[d]);
- if (rc < 0)
- goto error;
- rc = shrink_flow_z(flows_x[d], flows_y[d], NULL, NULL);
- if (rc < 0)
- goto error;
- rc = shrink_flow_z(flows_ix[d], flows_iy[d], NULL, NULL);
- if (rc < 0)
- goto error;
- }
- /* If any of the restricted pixles is 0 (marked),
- * the result must be 0 (marked) too.
- */
- quantize(marks[d], 0.9999, 1, 0);
- }
- if (d < scalexyz) {
- /* create x and y shrinked flow vectors for restriction and prolonging process */
- rc = shrink_xy(flows_x[d], flows_xy_x[d]);
- if (rc < 0)
- goto error;
- rc = shrink_xy(flows_y[d], flows_xy_y[d]);
- if (rc < 0)
- goto error;
- multiply(flows_xy_x[d], 0.5);
- multiply(flows_xy_y[d], 0.5);
- }
-
- /* generate values array */
- dimensions[0] = ww;
- dimensions[1] = hh;
- dimensions[2] = kk;
- values[d] = darrayCreate(3, dimensions);
- if (!values[d])
- goto error;
-
- /* generate array of neighbors */
-#ifndef TEST_NO_RENDER
- nb_list[d] = gen_neighbor(ww, hh, kk, flows_x[d], flows_y[d], flows_ix[d], flows_iy[d]);
- if (!nb_list[d])
- goto error;
- /* weighten */
- weighten(luminance, nb_list[d]);
-#endif
-
- last_kk = kk;
- if (d < scalexyz)
- kk = (kk+1)>>1;
- }
- darrayDestroy(luminance);
- luminance = NULL;
- printf("\n");
-
- dimensions[0] = w;
- dimensions[1] = h;
- dimensions[2] = k;
- init = darrayCreate(3, dimensions);
- if (!init)
- goto error;
-
- /* result array for all components */
- dimensions[0] = w;
- dimensions[1] = h;
- dimensions[2] = comp;
- dimensions[3] = k;
- result = darrayCreate(4, dimensions);
- if (!result)
- goto error;
-
-#ifdef TEST_PYRAMID
- dimensions[0] = w;
- dimensions[1] = h;
- dimensions[2] = k;
- luminance = darrayCreate(3, dimensions);
- for (z = 0; z < k; z++) {
- memcpy(darrayGetPr(luminance) + z*w*h, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
- memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
-#ifdef TEST_FLOW
- draw_flow(values[0], flows_x[0], flows_y[0], flows_ix[0], flows_iy[0]);
-#endif
-#ifdef TEST_MARKED
- paste(marks[0], values[0], 0, 0, 0);
-#endif
-
- }
- for (d = 1; d < max_depth; d++) {
- rc = shrink_xy(luminance, NULL);
- if (rc < 0)
- goto error;
- if (d-1 < scalexyz)
- shrink_z(luminance, NULL, flows_xy_x[d-1], flows_xy_y[d-1]);
- paste(luminance, values[0], 0, 0, 0);
-#ifdef TEST_FLOW
- draw_flow(values[0], flows_x[d], flows_y[d], flows_ix[d], flows_iy[d]);
-#endif
-#ifdef TEST_MARKED
- paste(marks[d], values[0], 0, 0, 0);
-#endif
- }
- darrayDestroy(luminance);
- luminance = NULL;
- /* copy test pyramid image as Y-component */
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-#else
- /* copy Y-component */
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
-#endif
-
- printf("-> apply component to grid and solve");
- for (c = 1; c < comp; c++) {
- int pixles, j;
- double *value_ptr, *init_ptr, *mark_ptr;
-
- printf(" #%d", c); fflush(stdout);
- /* apply component
- * - copy component into value array (level 0), otherwise 0
- * - use maked colors for init array, otherwise use 0
- */
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(image) + c*w*h + z*w*h*comp, w*h*sizeof(double));
- for (d = 1, ww=(w+1)>>1, hh=(h+1)>>1, kk=k; d < max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
- if (d-1 < scalexyz) /* scale kk, if last array was scaled */
- kk=(kk+1)>>1;
- value_ptr = darrayGetPr(values[d]);
- pixles = ww*hh*kk;
- for (j = 0; j < pixles; j++)
- *value_ptr++ = 0.0;
- }
- value_ptr = darrayGetPr(values[0]);
- init_ptr = darrayGetPr(init);
- mark_ptr = darrayGetPr(image_mark);
- pixles = w*h*k;
- for (j = 0; j < pixles; j++) {
- if (*mark_ptr++ > 0.5)
- *init_ptr++ = *value_ptr;
- else
- *init_ptr++ = 0.0;
- value_ptr++;
- }
-
-#ifndef TEST_NO_RENDER
- /* colorize component */
- int i;
- for (i = 0; i < outer_iter; i++) {
- rc = solve_mg(max_depth, inner_iter, 2, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
- if (rc < 0)
- goto error;
- rc = solve_mg((max_depth+1) / 2, inner_iter, 2, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
- if (rc < 0)
- goto error;
- rc = solve_mg(2, inner_iter, 2, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
- if (rc < 0)
- goto error;
- rc = solve_mg(1, inner_iter, 4, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
- if (rc < 0)
- goto error;
- }
-#endif
-
-#ifdef TEST_PYRAMID
- for (d = 1; d < max_depth; d++)
- paste(values[d], values[0], 0, 0, 0);
-#endif
-
- /* apply colorized component */
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
- }
- printf("\n");
-
-#ifdef TEST_SHRINK_XY
- for (c = 0; c < comp; c++) {
- darray_t *temp;
- int dimensions[3] = {w, h, k};
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
- dimensions[0] = (w+1)>>1;
- dimensions[1] = (h+1)>>1;
- temp = darrayCreate(3, dimensions);
- ASSERT(temp, "no memory");
- rc = shrink_xy(values[0], temp);
- ASSERT(rc >= 0, "no memory");
- paste(temp, values[0], 0, 0, 0);
- darrayDestroy(temp);
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
- }
-#endif
-
-#ifdef TEST_SHRINK_Z
- for (c = 0; c < comp; c++) {
- darray_t *temp;
- int dimensions[3] = {w, h, k};
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
- dimensions[2] = (k+1)>>1;
- temp = darrayCreate(3, dimensions);
- ASSERT(temp, "no memory");
- rc = shrink_z(values[0], temp, flows_x[0], flows_y[0]);
- ASSERT(rc >= 0, "no memory");
- paste(temp, values[0], 0, 0, 0);
- darrayDestroy(temp);
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
- }
-#endif
-
-#ifdef TEST_EXPAND_XY
- for (c = 0; c < comp; c++) {
- darray_t *temp;
- int dimensions[3] = {w, h, k};
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
- dimensions[0] = (w+1)>>1;
- dimensions[1] = (h+1)>>1;
- temp = darrayCreate(3, dimensions);
- ASSERT(temp, "no memory");
- rc = shrink_xy(values[0], temp);
- ASSERT(rc >= 0, "no memory");
- expand_xy(temp, values[0]);
- darrayDestroy(temp);
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
- }
-#endif
-
-#ifdef TEST_EXPAND_Z
- for (c = 0; c < comp; c++) {
- darray_t *temp;
- int dimensions[3] = {w, h, k};
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
- dimensions[2] = (k+1)>>1;
- temp = darrayCreate(3, dimensions);
- ASSERT(temp, "no memory");
- rc = shrink_z(values[0], temp, flows_x[0], flows_y[0]);
- ASSERT(rc >= 0, "no memory");
- rc = expand_z(temp, values[0], flows_x[0], flows_y[0]);
- ASSERT(rc >= 0, "no memory");
- darrayDestroy(temp);
- for (z = 0; z < k; z++)
- memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
- }
-#endif
-
- rc = 0;
-
-error:
- /* free memory */
- for (d = 0; d < max_depth; d++) {
- if (nb_list)
- free(nb_list[d]);
- if (marks)
- darrayDestroy(marks[d]);
- if (values)
- darrayDestroy(values[d]);
- if (flows_xy_x)
- darrayDestroy(flows_xy_x[d]);
- if (flows_xy_y)
- darrayDestroy(flows_xy_y[d]);
- if (flows_x)
- darrayDestroy(flows_x[d]);
- if (flows_y)
- darrayDestroy(flows_y[d]);
- if (flows_ix)
- darrayDestroy(flows_ix[d]);
- if (flows_iy)
- darrayDestroy(flows_iy[d]);
- }
- free(nb_list);
- free(marks);
- free(values);
- free(flows_xy_x);
- free(flows_xy_y);
- free(flows_x);
- free(flows_y);
- free(flows_ix);
- free(flows_iy);
- darrayDestroy(init);
- darrayDestroy(luminance);
- if (rc < 0)
- darrayDestroy(result);
-
- return (rc < 0) ? NULL : result;
-}
projects / colorize.git / blobdiff