/* Given two grey images, compare them using the algorithm described
in "An Image Signature for Any Kind of Image" by H. Chi Wong,
Marshall Bern and David Goldberg, 2002. */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdarg.h>
#include "similar-image.h"
#ifdef HEADER
#define SIZE 9
#define DIRECTIONS 8
/* A point in the grid. */
typedef struct point
{
double average_grey_level;
int d[DIRECTIONS];
}
point_t;
typedef int (*simage_error_channel_t) (void * s, const char * format, ...);
typedef struct simage
{
/* The width of the image in pixels. */
unsigned int width;
/* The height of the image in pixels. */
unsigned int height;
/* The image data. */
unsigned char * data;
/* The computed signature. */
char * signature;
/* The length of the signature. */
int signature_length;
/* The P-value for this image, see equation in article. */
unsigned int p;
/* The grid of values. */
point_t grid[SIZE*SIZE];
/* width / (SIZE + 1) */
double w10;
/* height / (SIZE + 1) */
double h10;
/* The number of times malloc/calloc were called related to this
object. */
int nmallocs;
simage_error_channel_t error_channel;
/* This contains a true value if we have actually loaded image
data, or a false value if not. This may be false if we just
loaded a signature rather than the image. */
unsigned int valid_image : 1;
/* The grid is already filled. */
unsigned int grid_filled : 1;
}
simage_t;
typedef enum
{
simage_ok,
/* malloc failed. */
simage_status_memory_failure,
/* x or y is outside the image dimensions. */
simage_status_bounds,
simage_status_bad_image,
/* Some upstream program did a stupid thing. */
simage_status_bad_logic,
/* */
simage_status_free_underflow,
/* */
simage_status_memory_leak,
}
simage_status_t;
typedef enum
{
much_darker = -2,
darker = -1,
same = 0,
lighter = 1,
much_lighter = 2,
}
comparison_t;
#endif /* def HEADER */
#define FAIL(test,status,message,...) { \
if (test) { \
if (s->error_channel) { \
(*s->error_channel) (s, "%s:%d: ", __FILE__, __LINE__); \
(*s->error_channel) (s, message, ## __VA_ARGS__); \
(*s->error_channel) (s, "\n"); \
} \
return simage_status_ ## status; \
} \
}
#define CALL(x) { \
simage_status_t status; \
status = x; \
if (status != simage_ok) { \
if (s->error_channel) { \
(*s->error_channel) (s, "%s:%d: ", __FILE__, __LINE__); \
(*s->error_channel) (s, "%s failed with status %d", \
#x, status); \
(*s->error_channel) (s, "\n"); \
} \
return status; \
} \
}
#define CHECK_XY(s,x,y) { \
FAIL (x > s->width || x < 0, bounds, \
"x coordinate %d is outside the image", x); \
FAIL (y > s->height || y < 0, bounds, \
"y coordinate %d is outside the image", y); \
}
/* Default place to print errors, if the user doesn't override this. */
static int
simage_default_error_channel (void * vs, const char * message, ...)
{
va_list va;
int chars;
va_start (va, message);
chars = vfprintf (stderr, message, va);
va_end (va);
return chars;
}
#define OUTSIDE -1
/* Given x and y coordinates, return the part of the grid which
corresponds to that. */
int x_y_to_entry (int x, int y)
{
int entry;
if (x < 0 || x >= SIZE) {
return OUTSIDE;
}
if (y < 0 || y >= SIZE) {
return OUTSIDE;
}
entry = y * SIZE + x;
if (entry < 0 || entry >= SIZE * SIZE) {
fprintf (stderr, "%s:%d: overflow %d\n", __FILE__, __LINE__, entry);
return OUTSIDE;
}
return entry;
}
simage_status_t
simage_dump (simage_t * s)
{
printf ("{\n");
printf ("\"width\":%d,\n", s->width);
printf ("\"height\":%d,\n", s->height);
printf ("\"p\":%d,\n", s->p);
printf ("\"dummy\":0\n");
printf ("}\n");
return simage_ok;
}
simage_status_t
simage_init (simage_t * s, unsigned int width, unsigned int height)
{
unsigned int p;
/* The minimum of the width and the height. */
unsigned int min_w_h;
s->data = calloc (width * height, sizeof (unsigned char));
CALL (simage_inc_nmallocs (s, s->data));
s->height = height;
s->width = width;
s->p = 2;
s->error_channel = & simage_default_error_channel;
min_w_h = width;
if (height < min_w_h) {
min_w_h = height;
}
p = (unsigned int) (floor (0.5 + ((double) min_w_h) / 20.0));
if (p > s->p) {
s->p = p;
}
// This contains a valid image data, although it is just black
// pixels at the moment.
s->valid_image = 1;
return simage_ok;
}
simage_status_t
simage_inc_nmallocs (simage_t * s, void * signature)
{
FAIL (! signature, memory_failure, "Out of memory");
s->nmallocs++;
return simage_ok;
}
simage_status_t
simage_dec_nmallocs (simage_t * s)
{
s->nmallocs--;
FAIL (s->nmallocs < 0, free_underflow,
"too many frees, %d should be 0.\n",
s->nmallocs);
return simage_ok;
}
/* Free all the memory associated with "s", except for "s" itself,
which is allocated by the user. */
simage_status_t
simage_free (simage_t * s)
{
if (s->data) {
free (s->data);
s->data = 0;
CALL (simage_dec_nmallocs (s));
}
if (s->signature) {
free (s->signature);
s->signature = 0;
CALL (simage_dec_nmallocs (s));
}
FAIL (s->nmallocs != 0, memory_leak,
"memory leak: %d should be 0.", s->nmallocs);
return simage_ok;
}
/* Set the pixel at "x", "y" to the value "grey". */
simage_status_t
simage_set_pixel (simage_t * s, int x, int y, unsigned char grey)
{
CHECK_XY (s, x, y);
s->data[y * s->width + x] = grey;
return simage_ok;
}
/* Compute the average intensity of the grid square at the coordinates
"i", "j" on the grid. This assumes that "s->w10" and "s->h10" have
already been computed. */
simage_status_t
simage_fill_entry (simage_t * s, int i, int j)
{
double total;
int px;
int py;
double xd;
double yd;
int x_min;
int y_min;
int x_max;
int y_max;
int size;
int entry;
int grey;
xd = (i + 1) * s->w10;
yd = (j + 1) * s->h10;
x_min = round (xd - s->p / 2.0);
y_min = round (yd - s->p / 2.0);
x_max = round (xd + s->p / 2.0);
y_max = round (yd + s->p / 2.0);
/* For very small images, these boundaries are sometimes
reached. */
if (y_max >= s->height) {
y_max = s->height - 1;
}
if (x_max >= s->width) {
x_max = s->width - 1;
}
if (x_min < 0) {
x_min = 0;
}
if (y_min < 0) {
y_min = 0;
}
total = 0.0;
for (py = y_min; py <= y_max; py++) {
FAIL (py < 0 || py >= s->height, bounds,
"overflow py=%d for i, j = (%d, %d)\n", py, i, j);
for (px = x_min; px <= x_max; px++) {
FAIL (px < 0 || px >= s->width, bounds,
"overflow px=%d for i, j = (%d, %d)\n", px, i, j);
total += s->data[py * s->width + px];
}
}
size = (x_max - x_min + 1) * (y_max - y_min + 1);
grey = (int) round (total / ((double) size));
FAIL (grey < 0 || grey >= 256, bounds,
"bad average grey value %d.", grey);
entry = x_y_to_entry (i, j);
FAIL (entry == OUTSIDE, bounds,
"bounds error with %d %d -> %d\n",
i, j, entry);
s->grid[entry].average_grey_level = grey;
return simage_ok;
}
/* Go around the image and make the average values for each of the
points on the grid. */
simage_status_t
simage_fill_entries (simage_t * s)
{
int i;
int j;
FAIL (s->width == 0 || s->height == 0, bad_image,
"empty image w/h %d/%d.\n",
s->width, s->height);
s->w10 = ((double) s->width) / ((double) (SIZE + 1));
s->h10 = ((double) s->height) / ((double) (SIZE + 1));
for (i = 0; i < SIZE; i++) {
for (j = 0; j < SIZE; j++) {
CALL (simage_fill_entry (s, i, j));
}
}
return simage_ok;
}
/* Given offsets xo and yo, return the array offset for the difference
array which corresponds to that. */
int xo_yo_to_direction (int xo, int yo)
{
int direction;
if (xo <= 0 && yo <= 0) {
direction = (xo + 1) + 3 * (yo + 1);
}
else {
// Adjust for not having a centre square, so +1, +1 is 7, not 8.
direction = (xo + 1) + 3 * (yo + 1) - 1;
}
return direction;
}
simage_status_t
direction_to_xo_yo (int direction, int * xo, int * yo)
{
if (direction < 3) {
* yo = -1;
* xo = direction - 1;
return simage_ok;
}
if (direction < 5) {
* yo = 0;
if (direction == 3) {
* xo = -1;
}
else if (direction == 4) {
* xo = 1;
}
else {
return simage_status_bounds;
}
return simage_ok;
}
if (direction < DIRECTIONS) {
* yo = 1;
* xo = direction - 6;
return simage_ok;
}
fprintf (stderr, "%s:%d: direction %d >= DIRECTIONS %d.\n",
__FILE__, __LINE__, direction, DIRECTIONS);
return simage_status_bounds;
}
int diff (int thisgrey, int thatgrey)
{
int d;
d = thisgrey - thatgrey;
if (d >= -2 && d <= 2) {
return same;
}
else if (d > 100) {
return much_darker;
}
else if (d > 2) {
return darker;
}
else if (d < -100) {
return much_lighter;
}
else if (d < -2) {
return lighter;
}
else {
fprintf (stderr, "%s:%d: mysterious d value %d\n",
__FILE__, __LINE__, d);
return same;
}
}
/* Make the difference between two adjoining points. */
simage_status_t
simage_make_point_diffs (simage_t * s, int x, int y)
{
int xo;
int yo;
int thisgrey;
int thisentry;
point_t * thispoint;
thisentry = x_y_to_entry (x, y);
/* Make 100% sure that we don't try to access outside the "grid" array
within "s". */
FAIL (thisentry == OUTSIDE, bounds, "entry outside grid %d %d %d\n",
x, y, thisentry);
thispoint = & s->grid[thisentry];
thisgrey = thispoint->average_grey_level;
for (xo = -1; xo <= 1; xo++) {
for (yo = -1; yo <= 1; yo++) {
int thatentry;
int direction;
int thatgrey;
if (xo == 0 && yo == 0) {
// Skip the middle square, since this would be the
// difference between us and ourselves.
continue;
}
thatentry = x_y_to_entry (x + xo, y + yo);
if (thatentry == OUTSIDE) {
// Skip entries which are outside the grid, which
// happens e.g. if x = 0 and xo = -1.
continue;
}
// Get the grey level of the other point
thatgrey = s->grid[thatentry].average_grey_level;
// turn xo, yo into an array offset "direction".
direction = xo_yo_to_direction (xo, yo);
// Put the difference into d[direction] of the current point.
thispoint->d[direction] = diff (thisgrey, thatgrey);
// fprintf (stderr, "# %d %d %d\n", thisentry, direction, thispoint->d[direction]);
}
}
return simage_ok;
}
simage_status_t
entry_to_x_y (int entry, int * x_ptr, int * y_ptr)
{
int x;
int y;
x = entry % SIZE;
y = entry / SIZE;
* x_ptr = x;
* y_ptr = y;
return simage_ok;
}
/* Make the array of differences between adjoining points. */
simage_status_t
simage_make_differences (simage_t * s)
{
int cell;
for (cell = 0; cell < SIZE * SIZE; cell++) {
int x;
int y;
CALL (entry_to_x_y (cell, & x, & y));
CALL (simage_make_point_diffs (s, x, y));
}
return simage_ok;
}
#define MAXDIM 10000
simage_status_t
simage_check_image (simage_t * s)
{
FAIL (s->width == 0 || s->height == 0, bad_image,
"empty image w/h %d/%d.\n",
s->width, s->height);
FAIL (s->width > MAXDIM || s->height > MAXDIM, bad_image,
"oversize image w/h %d/%d.\n",
s->width, s->height);
return simage_ok;
}
simage_status_t
simage_fill_grid (simage_t * s)
{
FAIL (s->grid_filled, bad_logic,
"double call to fill_grid.\n");
CALL (simage_check_image (s));
CALL (simage_fill_entries (s));
CALL (simage_make_differences (s));
s->grid_filled = 1;
return simage_ok;
}
simage_status_t
simage_diff (simage_t * s1, simage_t * s2, double * total_diff)
{
int total;
int total1;
int total2;
int cell;
total = 0;
total1 = 0;
total2 = 0;
for (cell = 0; cell < SIZE * SIZE; cell++) {
int direction;
for (direction = 0; direction < DIRECTIONS; direction++) {
int diff;
int s1cd;
int s2cd;
s1cd = s1->grid[cell].d[direction];
s2cd = s2->grid[cell].d[direction];
diff = s1cd - s2cd;
// Add the squares of the values to the totals.
total += diff * diff;
total1 += s1cd * s1cd;
total2 += s2cd * s2cd;
}
}
if (total1 == 0 && total2 == 0) {
*total_diff = 0.0;
return simage_ok;
}
*total_diff = ((double) total) / ((double)(total1 + total2));
return simage_ok;
}
/* Check whether this direction and cell point to another
cell or are outside the image. */
int inside (int cell, int direction)
{
int x;
int y;
int xo;
int yo;
int nextcell;
x = cell % SIZE;
y = cell / SIZE;
direction_to_xo_yo (direction, & xo, & yo);
nextcell = x_y_to_entry (x + xo, y + yo);
if (nextcell == OUTSIDE) {
return 0;
}
return 1;
}
simage_status_t
simage_allocate_signature (simage_t * s, int size)
{
s->signature = calloc (size + 1, sizeof (unsigned char));
CALL (simage_inc_nmallocs (s, s->signature));
return simage_ok;
}
/* Compute the signature of "s->data". If it has already been
computed, return without doing anything. */
simage_status_t
simage_signature (simage_t * s)
{
int cell;
int max_size;
if (s->signature) {
return simage_ok;
}
max_size = DIRECTIONS * SIZE * SIZE;
CALL (simage_allocate_signature (s, max_size));
s->signature_length = 0;
for (cell = 0; cell < SIZE * SIZE; cell++) {
int direction;
for (direction = 0; direction < DIRECTIONS; direction++) {
if (inside (cell, direction)) {
int value;
value = s->grid[cell].d[direction] + 2 + '0';
FAIL (value < '0' || value > '5', bounds,
"overflow %d at cell=%d direction=%d",
value, cell, direction);
s->signature[s->signature_length] = (char) value;
s->signature_length++;
FAIL (s->signature_length > max_size, bounds,
"signature length %d > max size %d",
s->signature_length, max_size);
}
}
}
return simage_ok;
}
simage_status_t
simage_fill_from_signature (simage_t * s, char * signature, int signature_length)
{
// Cell number
int c;
// Offset into signature.
int o;
CALL (simage_allocate_signature (s, signature_length));
s->signature_length = signature_length;
o = 0;
for (c = 0; c < SIZE * SIZE; c++) {
/* The direction. */
int d;
int x;
int y;
CALL (entry_to_x_y (c, & x, & y));
for (d = 0; d < DIRECTIONS; d++) {
if (inside (c, d)) {
//printf ("%d %d %d\n", c, d, o);
/* The value of this cell. */
int value;
FAIL (o >= signature_length, bounds,
"offset %d exceeded signature length %d.\n",
o, signature_length);
value = signature[o] - '0' - 2;
s->signature[o] = signature[o];
o++;
s->grid[c].d[d] = value;
}
}
}
s->signature[s->signature_length] = '\0';
return simage_ok;
}