for the merge with Jerome

Author: Pierre Paleo

comparisons/compare.c

@@ -3,31 +3,25 @@
  *
  * Usage : make && ./compare extrema_cpp.txt extrema_opencl.txt
  *
+ * Format of the keypoints in the files :
+ [  -3.41704249  368.90612793   14.99513245    2.78045988]
+ [  -6.37356377  449.18521118   69.4825592     2.53365469]
+  
  ********************************************************************
 */
 
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
+#include "compare.h"
 #define MAX_KP 300 //for a given octave 
 #define DIGITS 3 //for comparison precision (10e-DIGITS)
 
-typedef struct keypoint {
-	float p;
-	float r;
-	float c;
-	float s;
-} keypoint;
 
 /*
   Swap two keypoints pointers
 */
-void keypoint_swap(keypoint k1, keypoint k2) {
-	keypoint tmp_ptr = k1;
-	k1 = k2;
-	k2 = tmp_ptr; 
+void keypoint_swap(keypoint* k1, keypoint* k2) {
+	keypoint tmp_ptr = *k1;
+	*k1 = *k2;
+	*k2 = tmp_ptr; 
 }
 
 
@@ -55,7 +49,7 @@ int parse_keypoints(char* filename, keypoint* keypoints, unsigned int* total_key
 	FILE* stream = fopen(filename,"r+");
 	if (!stream) { printf("Error: Could not open file %s\n",filename); return -1; }
 	float p=0,r=0,c=0,s=0;
-	keypoint* kp = (keypoint*) calloc(1,sizeof(keypoint));
+	keypoint* kp = &(keypoints[0]);
 	unsigned int j= 1, k =0;
 	char str[511];
 	while (EOF != fscanf(stream,"%s",str)) {
@@ -73,17 +67,15 @@ int parse_keypoints(char* filename, keypoint* keypoints, unsigned int* total_key
 				case 0:
 					kp->s = atof(str);
 					keypoints[k] = *kp;
-					kp = (keypoint*) calloc(1,sizeof(keypoint));
+					kp = (keypoint*) &keypoints[k+1]; 
 					k++;
 					break;
 			} //end switch
 			j++;
 		} //end isdigit
 	} //end read loop
 	//puts a "end of keypoints" marker
-	kp = (keypoint*) calloc(1,sizeof(keypoint));
-	kp->r = -1.0;
-	keypoints[k] = *kp;
+	keypoints[k].r = -1.0;
 	*total_keypoints = k;
 	return 1;
 }
@@ -127,32 +119,71 @@ int main(int args, char* argv[]) {
 	}
 	printf("End of comparison -- %d/(%d,%d) keypoints matches\n",
 		kp_ok,total_keypoints_opencl,total_keypoints_cpp);
+	
+	keypoint* output = (keypoint*) calloc(total_keypoints_opencl-1,sizeof(keypoint));
+//	puts("before cut sort");
+	for(i=0;6 >= i; i++) printf("%f ",k_opencl[i].r);
+	//puts(""); puts("after cut sort");
+	merge_sort(k_opencl,0,6/*total_keypoints_opencl-1*/,output);
+	//for(i=0;30 > i; i++) printf("%f ",k_opencl[i].r);
+	puts("");
+	free(k_cpp);
+	free(k_opencl);
+	free(output);
 	return 1;
 
 
 }
 
-/*
-void cut_sort(keypoint* input, unsigned int start, unsigned int end) {
+//[start,...,end]
+void merge_sort(keypoint* input, unsigned int start, unsigned int end, keypoint* output) {
 	unsigned int len = end-start+1;
+	unsigned int middle = (end+start)/2;
 	if (len > 2) {
-		cut_sort(input, start, (end+start)/2);
-		cut_sort(input,(end+start)/2+1,end);
+		merge_sort(input, start, middle,output);
+		merge_sort(input,middle+1,end,output);
 	}
 	else {
-		if (len == 1) ; //return input[start];
+		if (len == 2) {
+			if (input[start].r > input[end].r) keypoint_swap(&input[start],&input[end]);
+		}
+	}
+/*	printf("Call to %d %d %d\n",start,middle,end);*/
+	if (len > 2) merge(input, output, start, middle, end);
+}
+/*
+	Merge 2 sorted lists
+*/
+void merge(keypoint* input, keypoint* output, unsigned int start, unsigned int middle, unsigned int end) {
+	int i1 = 0, i2 = 0, stop = 0;
+	printf("I am %d %d %d\n",start,middle,end);
+	while (!stop) {
+		if (input[start+i1].r > input[middle+1+i2].r) {
+			output[start+i1+i2] = input[middle+1+i2];
+			i2++;
+		}
 		else {
-			if (input[start].r > input[end].r) keypoint_swap(input[start],input[end]);
-			//else ;
+			output[start+i1+i2] = input[start+i1];
+			i1++;
 		}
-	}	
+		if (i1 == middle+1) stop = 1;
+		if (i2 == end+1) stop = 2;
+	}
+	int i;
+	if (stop == 1) //recopy the end of the 2nd list 
+		for (i = i2; end >= i; i++) output[start+i1+i] = input[middle+1+i];
+	if (stop == 2) //recopy the end of the 1st list
+		for (i = i1; middle >= i; i++) output[start+i2+i] = input[start+i];
 
 
+	puts("Merged :");
+	for (i = 0; i1+i2 >= i; i++) printf("%f ",output[start+i].r);
+	puts("");
+	
 }
 
 
-//TODO: print kp_opencl, see if kp_opencl[:].r have been sorted by parts of 2
-*/
+
 
 
 

compile

@@ -1,4 +1,4 @@
 export TD=dummy
 python setup.py install --install-lib ${TD}/libs --install-scripts ${TD}/scripts
 
-echo "N OUBLIE PAS DE INITIALISER PYTHONPATH VERS ${TD}/libs"
+echo "N OUBLIE PAS DE INITIALISER PYTHONPATH VERS -path absolu- de ${TD}/libs"

openCL/image.cl

@@ -385,6 +385,7 @@ __kernel void interp_keypoint(
  * @param grad: Pointer to global memory with gradient norm previously calculated
  * @param ori: Pointer to global memory with gradient orientation previously calculated
  * @param counter: Pointer to global memory with actual number of keypoints previously found
+ * @param hist: Pointer to shared memory with histogram (36 values per thread)
  * @param octsize: initially 1 then twiced at each octave
  * @param OriSigma : a SIFT parameter, default is 1.5. Warning : it is not "InitSigma".
  * @param nb_keypoints : maximum number of keypoints
@@ -400,13 +401,28 @@ par.OriHistThresh = 0.8;
 -replace "36" by an external paramater ?
 -replace "0.8" by an external parameter ?
 
+TODO:
+-Memory optimization
+	--Use less registers (re-use, calculation instead of assignation)
+	--Use local memory for float histogram[36]
+-Speed-up
+	--Less access to global memory (k.s1 is OK because this is a register)
+	--leave the loops as soon as possible
+	--Avoid divisions
+
+TODO: QUESTION :
+keypoint k = keypoints[gid0];
+keypoint is float4.
+Is this a copy (register) or a pointer ? would explain bad performances...
+
 */
 
 __kernel void orientation_assignment(
 	__global keypoint* keypoints,
 	__global float* grad, 
 	__global float* ori,
 	__global int* counter,
+	__local float* hist,
 	int octsize,
 	float OriSigma, //WARNING: (1.5), it is not "InitSigma (=1.6)"
 	int nb_keypoints,
@@ -420,33 +436,37 @@ __kernel void orientation_assignment(
 	if (keypoints_start <= gid0 && gid0 < keypoints_end) { //do not use *counter, for it will be modified below
 		keypoint k = keypoints[gid0];
 		if (k.s1 != -1.0f) { //if the keypoint is valid 
+			
+			//Local memory memset
+			for (int i=0; i < 36; i++)
+				hist[36*gid0+i] = 0.0f;
+				
 			int	bin, prev, next;
-			int old;
+			int old; //counter value
 			float distsq, gval, angle, interp=0.0;
-			float hist[36] = { 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f, 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f, 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f};
-			int	row = (int) (k.s1 + 0.5),
-				col = (int) (k.s2 + 0.5);
+			//int	row = (int) (k.s1 + 0.5),
+			//	col = (int) (k.s2 + 0.5);
 
 			/* Look at pixels within 3 sigma around the point and sum their
 			  Gaussian weighted gradient magnitudes into the histogram. */
 
 			float sigma = OriSigma * k.s3;
 			int	radius = (int) (sigma * 3.0);
-			int rmin = MAX(0,row - radius);
-			int cmin = MAX(0,col - radius);
-			int rmax = MIN(row + radius,grad_height - 2);
-			int cmax = MIN(col + radius,grad_width - 2);
+			int rmin = MAX(0,((int) (k.s1 + 0.5)) - radius);
+			int cmin = MAX(0,((int) (k.s2 + 0.5)) - radius);
+			int rmax = MIN(((int) (k.s1 + 0.5)) + radius,grad_height - 2);
+			int cmax = MIN(((int) (k.s2 + 0.5)) + radius,grad_width - 2);
 			int i,j,r,c;
 			for (r = rmin; r <= rmax; r++) {
 				for (c = cmin; c <= cmax; c++) {
 
 					gval = grad[r*grad_width+c];
-					distsq = (r-k.s1)*(r-k.s1) + (c-k.s2)*(c-k.s2);
+					distsq = (r-k.s1,2)*(r-k.s1,2) + (c-k.s2,2)*(c-k.s2,2);
 
 					if (gval > 0.0f  &&  distsq < ((float) (radius*radius)) + 0.5f) {
 						/* Ori is in range of -PI to PI. */
 						angle = ori[r*grad_width+c];
-						bin = (int) (36 * (angle + M_PI_F + 0.001f) / (2.0f * M_PI_F)); //FIXME: why this offset ?
+						bin = (int) (36 * (angle + M_PI_F + 0.001f) / (2.0f * M_PI_F)); //why this offset ?
 						if (bin >= 0 && bin <= 36) {
 							bin = MIN(bin, 35);
 							hist[bin] += exp(- distsq / (2.0f*sigma*sigma)) * gval;
@@ -485,22 +505,23 @@ __kernel void orientation_assignment(
 				hist[next] = -hist[next];
 			}
 			interp = 0.5f * (hist[prev] - hist[next]) / (hist[prev] - 2.0f * maxval + hist[next]);
-			angle = 2.0f * M_PI_F * (argmax + 0.5f + interp) / 36 - M_PI_F;
-		
+			angle = 2.0f * M_PI_F * (argmax + 0.5f + interp) * 0.027778f - M_PI_F; //1/36 = 0.027777777777777776
 
+			/*
+				Re-arrange coordinates to be coherent with sift.cpp
+			*/
 			k.s0 = k.s2; //c
 			k.s1 = k.s1; //r
 			k.s2 = k.s3; //sigma
 			k.s3 = angle;		   //angle
-		
 			keypoints[gid0] = k;
 
 			/*
 				An orientation is now assigned to our current keypoint.
 				We can create new keypoints of same (x,y,sigma) but a different angle.
 			 	For every local peak in histogram, every peak of value >= 80% of maxval generates a new keypoint	
 			*/
-		
+			//TODO: use k instead of k2 for memory ?
 			keypoint k2 = 0.0; k2.s0 = k.s0; k2.s1 = k.s1; k2.s2 = k.s2;
 			for (i = 0; i < 36; i++) {
 				prev = (i == 0 ? 36 -1 : i - 1);

test/test_image.py

@@ -269,10 +269,13 @@ def test_orientation(self):
         keypoints_start = numpy.int32(0)
         keypoints_end = numpy.int32(actual_nb_keypoints)
         counter = pyopencl.array.to_device(queue, keypoints_end) #actual_nb_keypoints)
+        #shared memory : 36-bin histograms
+        local_size = (keypoints_end-keypoints_start+1)*36*4
+        local_mem = pyopencl.LocalMemory(local_size)
 
         t0 = time.time()
         k1 = self.program.orientation_assignment(queue, shape, wg, 
-        	gpu_keypoints.data, gpu_grad.data, gpu_ori.data, counter.data,
+        	gpu_keypoints.data, gpu_grad.data, gpu_ori.data, counter.data, local_mem,
         	octsize, orisigma, nb_keypoints, keypoints_start, keypoints_end, grad_width, grad_height)    	
         res = gpu_keypoints.get()
         cnt = counter.get()