Re[2]: [XviD-devel] motion estimation for B/P/I decision

[Marco Al]

This is a multi-part message in MIME format.

------=_NextPart_000_0007_01C260E0.1FEAAF50
Content-Type: text/plain;
	charset="iso-8859-1"
Content-Transfer-Encoding: 7bit

Marc FD wrote:

> @MarcoAI
> do you have any C implementation of Zhou Wang's structural distortion
> measure ??

Here is a working one (ignore the C++ stuff around it, thats just for testing it
with the 512*512 GIFs from Zhou Whang's site).

Marco

------=_NextPart_000_0007_01C260E0.1FEAAF50
Content-Type: text/plain;
	name="distort.cpp"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
	filename="distort.cpp"

// distort.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "atlimage.h"
#include <iostream.h>

#define BLOCK_SIZE 8
#define N BLOCK_SIZE * BLOCK_SIZE

float
image_structural_distortion(unsigned char *img1,
		   unsigned char *img2,
		   int stride,
		   int width,
		   int height)
{
	int x, y, img1_sum_temp, img2_sum_temp, img1_sq_sum_temp, =
img2_sq_sum_temp,
		img12_sum_temp, img1_sq_sum, img2_sq_sum, img12_sum, img12_sum_mul,
		img12_sum_sq_sum;
	short img1_sum, img2_sum;
	float numerator, denominator1, denominator, q, structural_distortion =
=3D 0;

	//alignment wont do a hell of a lot of good given the way we access =
this ...
	short *img1_sum_ptr, *img1_sum_acc =3D (short *)malloc((height + 1) * =
(width + 1) * 2);
	short *img2_sum_ptr, *img2_sum_acc =3D (short *)malloc((height + 1) * =
(width + 1) * 2);
	int *img1_sq_sum_ptr, *img1_sq_sum_acc =3D (int *)malloc((height + 1) * =
(width + 1) * 4);
	int *img2_sq_sum_ptr, *img2_sq_sum_acc =3D (int *)malloc((height + 1) * =
(width + 1) * 4);
	int *img12_sum_ptr, *img12_sum_acc =3D (int *)malloc((height + 1) * =
(width + 1) * 4);

	img1_sum_ptr =3D img1_sum_acc;
	img2_sum_ptr =3D img2_sum_acc;
	img1_sq_sum_ptr =3D img1_sq_sum_acc;
	img2_sq_sum_ptr =3D img2_sq_sum_acc;
	img12_sum_ptr =3D img12_sum_acc;

	for (x =3D 0; x < width + 1; x++) {
		*(img1_sum_acc + x) =3D 0;
		*(img2_sum_acc + x) =3D 0;
		*(img1_sq_sum_acc + x) =3D 0;
		*(img2_sq_sum_acc + x) =3D 0;
		*(img12_sum_acc + x) =3D 0;
	}

	for (y =3D 1; y < height + 1; y++) {
		*(img1_sum_acc + y * (width + 1)) =3D 0;
		*(img2_sum_acc + y * (width + 1)) =3D 0;
		*(img1_sq_sum_acc + y * (width + 1)) =3D 0;
		*(img2_sq_sum_acc + y * (width + 1)) =3D 0;
		*(img12_sum_acc + y * (width + 1)) =3D 0;
	}

	img1_sum_acc +=3D width + 2;
	img2_sum_acc +=3D width + 2;
	img1_sq_sum_acc +=3D width + 2;
	img2_sq_sum_acc +=3D width + 2;
	img12_sum_acc +=3D width + 2;

	for (y =3D 0; y < height; y++) {
		img1_sum_temp =3D 0;
		img2_sum_temp =3D 0;
		img1_sq_sum_temp =3D 0;
		img2_sq_sum_temp =3D 0;
		img12_sum_temp =3D 0;

		for (x =3D 0; x < width; x++) {
			img1_sum_temp +=3D  *(img1 + x);
			*(img1_sum_acc + x) =3D *(img1_sum_acc + x - (width + 1)) + =
img1_sum_temp;
			img2_sum_temp +=3D *(img2 + x);
			*(img2_sum_acc + x) =3D *(img2_sum_acc + x - (width + 1)) + =
img2_sum_temp;
			img1_sq_sum_temp +=3D *(img1 + x) * *(img1 + x);
			*(img1_sq_sum_acc + x) =3D *(img1_sq_sum_acc + x - (width + 1)) + =
img1_sq_sum_temp;
			img2_sq_sum_temp +=3D *(img2 + x) * *(img2 + x);
			*(img2_sq_sum_acc + x) =3D *(img2_sq_sum_acc + x - (width + 1)) + =
img2_sq_sum_temp;
			img12_sum_temp +=3D *(img1 + x) * *(img2 + x);
			*(img12_sum_acc + x) =3D *(img12_sum_acc + x - (width + 1)) + =
img12_sum_temp;
		}

		img1 +=3D stride;
		img2 +=3D stride;
		img1_sum_acc +=3D width + 1;
		img2_sum_acc +=3D width + 1;
		img1_sq_sum_acc +=3D width + 1;
		img2_sq_sum_acc +=3D width + 1;
		img12_sum_acc +=3D width + 1;
	}

	img1_sum_acc =3D img1_sum_ptr;
	img2_sum_acc =3D img2_sum_ptr;
    img1_sq_sum_acc =3D img1_sq_sum_ptr;
	img2_sq_sum_acc =3D img2_sq_sum_ptr;
    img12_sum_acc =3D img12_sum_ptr;

	for (y =3D 0; y < height - BLOCK_SIZE + 1; y++) {
		for (x =3D 0; x < width - BLOCK_SIZE + 1; x++) {
			img1_sum =3D *(img1_sum_acc + x + BLOCK_SIZE + BLOCK_SIZE * (width + =
1)) -
                *(img1_sum_acc + x + BLOCK_SIZE) -
				*(img1_sum_acc + x + BLOCK_SIZE * (width + 1)) +
				*(img1_sum_acc + x);
			img2_sum =3D *(img2_sum_acc + x + BLOCK_SIZE + BLOCK_SIZE * (width + =
1)) -
				*(img2_sum_acc + x + BLOCK_SIZE) -
				*(img2_sum_acc + x + BLOCK_SIZE * (width + 1)) +
				*(img2_sum_acc + x);
			img1_sq_sum =3D *(img1_sq_sum_acc + x + BLOCK_SIZE + BLOCK_SIZE * =
(width + 1)) -
				*(img1_sq_sum_acc + x + BLOCK_SIZE) -
				*(img1_sq_sum_acc + x + BLOCK_SIZE * (width + 1)) +
				*(img1_sq_sum_acc + x);
			img2_sq_sum =3D *(img2_sq_sum_acc + x + BLOCK_SIZE + BLOCK_SIZE * =
(width + 1)) -
				*(img2_sq_sum_acc + x + BLOCK_SIZE) -
				*(img2_sq_sum_acc + x + BLOCK_SIZE * (width + 1)) +
				*(img2_sq_sum_acc + x);
			img12_sum =3D *(img12_sum_acc + x + BLOCK_SIZE + BLOCK_SIZE * (width =
+ 1)) -
				*(img12_sum_acc + x + BLOCK_SIZE) -
				*(img12_sum_acc + x + BLOCK_SIZE * (width + 1)) +
				*(img12_sum_acc + x);
			img12_sum_mul =3D img1_sum * img2_sum;
			/* the matlab equivalent for the next name was silly, yes even =
sillier than these ones, so instead
			of sq sum mul I used sum sq sum */
			img12_sum_sq_sum =3D img1_sum * img1_sum + img2_sum * img2_sum;

			numerator =3D (float)4 * ((float)N * (float)img12_sum - =
(float)img12_sum_mul) * (float)img12_sum_mul;
			denominator1 =3D (float)N * ((float)img1_sq_sum + (float)img2_sq_sum) =
- (float)img12_sum_sq_sum;
			denominator =3D denominator1 * (float)img12_sum_sq_sum;

			q =3D 1;
			if ((denominator1 =3D=3D 0) && (img12_sum_sq_sum !=3D 0))
				q =3D 2 * (float)img12_sum_mul / (float)img12_sum_sq_sum;
			if (denominator !=3D 0)
				q =3D numerator / denominator;
		=09
			structural_distortion +=3D q;
		}

		img1_sum_acc +=3D width + 1;
		img2_sum_acc +=3D width + 1;
		img1_sq_sum_acc +=3D width + 1;
		img2_sq_sum_acc +=3D width + 1;
		img12_sum_acc +=3D width + 1;
	}
=09
	structural_distortion /=3D (width - BLOCK_SIZE + 1) * (height - =
BLOCK_SIZE + 1);
	return structural_distortion;
}

int _tmain(int argc, _TCHAR* argv[])
{
    CImage im1, im2;
    im1.Load(argv[1]);=20
    im2.Load(argv[2]);

	unsigned char *img1 =3D (unsigned char *)im1.GetBits();
	unsigned char *img2 =3D (unsigned char *)im2.GetBits();

	img1 -=3D 511*512;
	img2 -=3D 511*512;

	cout << "Q: " << image_structural_distortion(img1, img2, 512, 512, =
512);

	return 0;
}

------=_NextPart_000_0007_01C260E0.1FEAAF50--