2 * @file interferometer.c
3 * @purpose Implementation of interferometer related functions
8 #include "interferometer.h"
11 /** Buffer for storing image data. Stored as a single intensity value for the laser light **/
12 static CvMat * g_data = NULL;
15 /** Camera capture pointer **/
16 static CvCapture * g_capture = NULL;
23 // For testing purposes
24 double test_omega = 0.05;
25 double test_angle = PI/2;
26 double test_noise[] = {0,0,0.02};
27 double test_phase = 0;
28 double test_intensity = 1.0;
31 static void Interferometer_TestSinusoid()
33 if (g_capture == NULL)
35 g_capture = cvCreateCameraCapture(0);
38 // Get image from camera
39 IplImage * img = cvQueryFrame(g_capture);
43 CvMat * background = cvGetMat(img, &stub, 0, 0);
44 // ... Honestly, I have no idea what the "stub" is for
48 g_data = cvCreateMat(background->rows, background->cols, CV_32FC3);
51 //cvShowImage("background", background);
53 for (int x = 0; x < g_data->cols-1; ++x)
55 for (int y = 0; y < g_data->rows-1; ++y)
57 // Calculate pure sine in test direction
58 double r = x*cos(test_angle) + y*sin(test_angle);
59 double value = 0.5*test_intensity*(1+sin(test_omega*r + test_phase));
62 s.val[0] = 0; s.val[1] = 0; s.val[2] = value;
63 CvScalar b = cvGet2D(background, y, x);
66 // Add noise & background image
68 // Get the order the right way round
72 for (int i = 0; i < 3; ++i)
75 s.val[i] += (rand() % 1000) * 1e-3 * test_noise[i];
76 s.val[i] += b.val[i] / 255; // Camera image is 0-255
79 //printf("set %d,%d\n", x, y);
83 cvSet2D(g_data,y, x, s);
92 * Get an image from the Interferometer
94 static void Interferometer_GetImage()
98 Interferometer_TestSinusoid();
99 //TODO: Implement camera
104 * Initialise the Interferometer
106 void Interferometer_Init()
109 // Make an initial reading (will allocate memory the first time only).
110 Interferometer_Read(1);
114 * Cleanup Interferometer stuff
116 void Interferometer_Cleanup()
119 cvReleaseMat(&g_data);
121 if (g_capture != NULL)
122 cvReleaseCapture(&g_capture);
127 * Read the interferometer; gets the latest image, processes it, spits out a single number
128 * @param samples - Number of columns to scan (increasing will slow down performance!)
129 * @returns Value proportional to the change in interferometer path length since the last call to this function
131 double Interferometer_Read(int samples)
137 // Get the latest image
138 Interferometer_GetImage();
139 // Frequency of the sinusoid
140 static double omega = 0;
141 // Stores locations of nodes
142 static int nodes[MAXNODES];
144 static double phase = 0;
147 double cur_phase = 0;
149 int xstep = g_data->cols / (samples+1);
151 // Used for testing to see where the nodes are identified
152 // (Can't modify g_data)
153 static CvMat * test_overlay = NULL;
154 if (test_overlay == NULL)
156 // Creates a memory leak; don't do this in the final version!
157 test_overlay = cvCreateMat(g_data->rows, g_data->cols, CV_32FC3);
159 cvZero(test_overlay);
160 //cvCopy(g_data, test_overlay, NULL);
162 // For each column to sample
163 for (int x = xstep; x < g_data->cols; x += xstep)
166 double avg = 0.5; //TODO: Calculate from image
167 double threshold_dif = 0; //TODO: Pick this value
172 for (int y = 1; y < g_data->rows-2 && num_nodes < MAXNODES; ++y)
174 if (num_nodes == 0 || abs(nodes[num_nodes-1] - y) > 1)
176 // A "node" is defined where the ajacent points are on opposite sides of the avg
177 double ldif = INTENSITY(cvGet2D(g_data, y-1,x)) - avg;
178 double rdif = INTENSITY(cvGet2D(g_data, y+1,x)) - avg;
180 // If that is the case, the product of the differences will be negative
181 if (ldif * rdif < -threshold_dif)
184 nodes[num_nodes++] = y;
186 // Put a white line on the overlay to indicate the node was found
187 for (int xx = 0; xx < g_data->cols; ++xx)
189 CvScalar s; // = cvGet2D(g_data, y, xx);
190 s.val[0] = 1; s.val[1] = 1; s.val[2] = 1;
191 cvSet2D(test_overlay, y, xx, s);
197 // Insufficient nodes found to continue
204 // Estimate angular frequency from two nodes TODO: Average between nodes
205 double slice_omega = (PI*(num_nodes-1)) / (nodes[num_nodes-1] - nodes[0]);
206 //printf("SLICE: %f vs %f\n", slice_omega, test_omega);
208 double slice_phase = 0;
209 for (int i = 0; i < num_nodes; ++i)
211 double this_phase = ((double)(i)*PI - slice_omega * nodes[i]);
212 //printf("Node %d gives phase %f\n", i, this_phase);
213 slice_phase += this_phase;
215 slice_phase /= num_nodes;
216 cur_phase += slice_phase;
219 // Average over samples
223 cur_phase /= samples;
227 // Get phase change since last call, save current phase
228 double result = (cur_phase - phase);
231 if (abs(result) > 0.5*PI)
240 // Display the image with lines to indicate results of data processing
241 //cvShowImage("nodes", test_overlay);
243 cvAdd(test_overlay, g_data, test_overlay, NULL);
244 cvShowImage("overlay", test_overlay);
251 * For testing purposes
253 int main(int argc, char ** argv)
255 //cvNamedWindow( "display", CV_WINDOW_AUTOSIZE );// Create a window for display.
256 gettimeofday(&start, NULL);
258 //Interferometer_Read(1);
259 //exit(EXIT_SUCCESS);
261 Interferometer_Init();
264 double last_phase = test_phase;
271 gettimeofday(&now, NULL);
272 time = TIMEVAL_DIFF(now, start);
274 test_phase = 0.5*PI*sin(time);
277 double delta = Interferometer_Read(1);
280 printf("%f\t%f\t%f\t%f\t%f\n", time, test_phase - last_phase, test_phase, delta, sum);
282 last_phase = test_phase;