view release on metacpan or  search on metacpan

ccv-src/lib/ccv_sift.c  view on Meta::CPAN

	if (x >= 0 || (float)xi == x)
		return xi;
	return xi - 1;
}

#define EXPN_SZ  256 /* fast_expn table size */
#define EXPN_MAX 25.0 /* fast_expn table max */
static double _ccv_expn_tab[EXPN_SZ + 1]; /* fast_expn table */
static int _ccv_expn_init = 0;

static inline double _ccv_expn(double x)
{
	double a, b, r;
	int i;
	assert(0 <= x && x <= EXPN_MAX);
	if (x > EXPN_MAX)
		return 0.0;
	x *= EXPN_SZ / EXPN_MAX;
	i = (int)x;
	r = x - i;
	a = _ccv_expn_tab[i];
	b = _ccv_expn_tab[i + 1];
	return a + r * (b - a);
}

static void _ccv_precomputed_expn()
{
	int i;
	for(i = 0; i < EXPN_SZ + 1; i++)
		_ccv_expn_tab[i] = exp(-(double)i * (EXPN_MAX / EXPN_SZ));
	_ccv_expn_init = 1;
}

void ccv_sift(ccv_dense_matrix_t* a, ccv_array_t** _keypoints, ccv_dense_matrix_t** _desc, int type, ccv_sift_param_t params)
{
	assert(CCV_GET_CHANNEL(a->type) == CCV_C1);
	ccv_dense_matrix_t** g = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels + 1) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	memset(g, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels + 1) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	ccv_dense_matrix_t** dog = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels - 1) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	memset(dog, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels - 1) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	ccv_dense_matrix_t** th = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	memset(th, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	ccv_dense_matrix_t** md = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	memset(md, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x ? params.noctaves + 1 : params.noctaves));
	if (params.up2x)
	{
		g += params.nlevels + 1;
		dog += params.nlevels - 1;
		th += params.nlevels - 3;
		md += params.nlevels - 3;
	}
	ccv_array_t* keypoints = *_keypoints;
	int custom_keypoints = 0;
	if (keypoints == 0)
		keypoints = *_keypoints = ccv_array_new(sizeof(ccv_keypoint_t), 10, 0);
	else
		custom_keypoints = 1;
	int i, j, k, x, y;
	double sigma0 = 1.6;
	double sigmak = pow(2.0, 1.0 / (params.nlevels - 3));
	double dsigma0 = sigma0 * sigmak * sqrt(1.0 - 1.0 / (sigmak * sigmak));
	if (params.up2x)
	{
		ccv_sample_up(a, &g[-(params.nlevels + 1)], 0, 0, 0);
		/* since there is a gaussian filter in sample_up function already,
		 * the default sigma for upsampled image is sqrt(2) */
		double sd = sqrt(sigma0 * sigma0 - 2.0);
		ccv_blur(g[-(params.nlevels + 1)], &g[-(params.nlevels + 1) + 1], CCV_32F | CCV_C1, sd);
		ccv_matrix_free(g[-(params.nlevels + 1)]);
		for (j = 1; j < params.nlevels; j++)
		{
			sd = dsigma0 * pow(sigmak, j - 1);
			ccv_blur(g[-(params.nlevels + 1) + j], &g[-(params.nlevels + 1) + j + 1], 0, sd);
			ccv_subtract(g[-(params.nlevels + 1) + j + 1], g[-(params.nlevels + 1) + j], (ccv_matrix_t**)&dog[-(params.nlevels - 1) + j - 1], 0);
			if (j > 1 && j < params.nlevels - 1)
				ccv_gradient(g[-(params.nlevels + 1) + j], &th[-(params.nlevels - 3) + j - 2], 0, &md[-(params.nlevels - 3) + j - 2], 0, 1, 1);
			ccv_matrix_free(g[-(params.nlevels + 1) + j]);
		}
		ccv_matrix_free(g[-1]);
	}
	double sd = sqrt(sigma0 * sigma0 - 0.25);
	g[0] = a;
	/* generate gaussian pyramid (g, dog) & gradient pyramid (th, md) */
	ccv_blur(g[0], &g[1], CCV_32F | CCV_C1, sd);
	for (j = 1; j < params.nlevels; j++)
	{
		sd = dsigma0 * pow(sigmak, j - 1);
		ccv_blur(g[j], &g[j + 1], 0, sd);
		ccv_subtract(g[j + 1], g[j], (ccv_matrix_t**)&dog[j - 1], 0);
		if (j > 1 && j < params.nlevels - 1)
			ccv_gradient(g[j], &th[j - 2], 0, &md[j - 2], 0, 1, 1);
		ccv_matrix_free(g[j]);
	}
	ccv_matrix_free(g[params.nlevels]);
	for (i = 1; i < params.noctaves; i++)
	{
		ccv_sample_down(g[(i - 1) * (params.nlevels + 1)], &g[i * (params.nlevels + 1)], 0, 0, 0);
		if (i - 1 > 0)
			ccv_matrix_free(g[(i - 1) * (params.nlevels + 1)]);
		sd = sqrt(sigma0 * sigma0 - 0.25);
		ccv_blur(g[i * (params.nlevels + 1)], &g[i * (params.nlevels + 1) + 1], CCV_32F | CCV_C1, sd);
		for (j = 1; j < params.nlevels; j++)
		{
			sd = dsigma0 * pow(sigmak, j - 1);
			ccv_blur(g[i * (params.nlevels + 1) + j], &g[i * (params.nlevels + 1) + j + 1], 0, sd);
			ccv_subtract(g[i * (params.nlevels + 1) + j + 1], g[i * (params.nlevels + 1) + j], (ccv_matrix_t**)&dog[i * (params.nlevels - 1) + j - 1], 0);
			if (j > 1 && j < params.nlevels - 1)
				ccv_gradient(g[i * (params.nlevels + 1) + j], &th[i * (params.nlevels - 3) + j - 2], 0, &md[i * (params.nlevels - 3) + j - 2], 0, 1, 1);
			ccv_matrix_free(g[i * (params.nlevels + 1) + j]);
		}
		ccv_matrix_free(g[i * (params.nlevels + 1) + params.nlevels]);
	}
	ccv_matrix_free(g[(params.noctaves - 1) * (params.nlevels + 1)]);
	if (!custom_keypoints)
	{
		/* detect keypoint */
		for (i = (params.up2x ? -1 : 0); i < params.noctaves; i++)
		{
			double s = pow(2.0, i);
			int rows = dog[i * (params.nlevels - 1)]->rows;
			int cols = dog[i * (params.nlevels - 1)]->cols;
			for (j = 1; j < params.nlevels - 2; j++)
			{
				float* bf = dog[i * (params.nlevels - 1) + j - 1]->data.f32 + cols;
				float* cf = dog[i * (params.nlevels - 1) + j]->data.f32 + cols;
				float* uf = dog[i * (params.nlevels - 1) + j + 1]->data.f32 + cols;
				for (y = 1; y < rows - 1; y++)
				{
					for (x = 1; x < cols - 1; x++)
					{
						float v = cf[x];
#define locality_if(CMP, SGN) \
	(v CMP ## = SGN params.peak_threshold && v CMP cf[x - 1] && v CMP cf[x + 1] && \
	 v CMP cf[x - cols - 1] && v CMP cf[x - cols] && v CMP cf[x - cols + 1] && \
	 v CMP cf[x + cols - 1] && v CMP cf[x + cols] && v CMP cf[x + cols + 1] && \
	 v CMP bf[x - 1] && v CMP bf[x] && v CMP bf[x + 1] && \
	 v CMP bf[x - cols - 1] && v CMP bf[x - cols] && v CMP bf[x - cols + 1] && \
	 v CMP bf[x + cols - 1] && v CMP bf[x + cols] && v CMP bf[x + cols + 1] && \
	 v CMP uf[x - 1] && v CMP uf[x] && v CMP uf[x + 1] && \
	 v CMP uf[x - cols - 1] && v CMP uf[x - cols] && v CMP uf[x - cols + 1] && \
	 v CMP uf[x + cols - 1] && v CMP uf[x + cols] && v CMP uf[x + cols + 1])
						if (locality_if(<, -) || locality_if(>, +))
						{
							ccv_keypoint_t kp;
							int ix = x, iy = y;
							double score = -1;
							int cvg = 0;
							int offset = ix + (iy - y) * cols;
							/* iteratively converge to meet subpixel accuracy */
							for (k = 0; k < 5; k++)
							{
								offset = ix + (iy - y) * cols;
								float N9[3][9] = { { bf[offset - cols - 1], bf[offset - cols], bf[offset - cols + 1],
													 bf[offset - 1], bf[offset], bf[offset + 1],
													 bf[offset + cols - 1], bf[offset + cols], bf[offset + cols + 1] },
												   { cf[offset - cols - 1], cf[offset - cols], cf[offset - cols + 1],
													 cf[offset - 1], cf[offset], cf[offset + 1],
													 cf[offset + cols - 1], cf[offset + cols], cf[offset + cols + 1] },
												   { uf[offset - cols - 1], uf[offset - cols], uf[offset - cols + 1],
													 uf[offset - 1], uf[offset], uf[offset + 1],
													 uf[offset + cols - 1], uf[offset + cols], uf[offset + cols + 1] } };
								score = _ccv_keypoint_interpolate(N9, ix, iy, j, &kp);
								if (kp.x >= 1 && kp.x <= cols - 2 && kp.y >= 1 && kp.y <= rows - 2)
								{