1 /**********
2 *
3 * This file is part of the TINA Open Source Image Analysis Environment
4 * henceforth known as TINA
5 *
6 * TINA is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as
8 * published by the Free Software Foundation.
9 *
10 * TINA is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with TINA; if not, write to the Free Software Foundation, Inc.,
17 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 **********
20 *
21 * Program : TINA
22 * File : $Source: /home/tina/cvs/tina-libs/tina/vision/visCorr_correlate.c,v $
23 * Date : $Date: 2008/10/30 09:43:59 $
24 * Version : $Revision: 1.2 $
25 * CVS Id : $Id: visCorr_correlate.c,v 1.2 2008/10/30 09:43:59 neil Exp $
26 *
27 * Author : S. Crossley, R. Lane (see thesis Univ. Sheffield U.K.), NAT.
28 */
29 /**
30 * @file stereo correlation algorithms
31 * @brief Stretch correlation algorithms (see Richard Lanes's Thesis Univ. Sheffiled).
32 *
33 *********
34 */
35
36 #include "visCorr_correlate.h"
37
38 #include <math.h>
39 #include <float.h>
40 #include <limits.h>
41 #include <tina/sys/sysDef.h>
42 #include <tina/sys/sysPro.h>
43 #include <tina/math/mathDef.h>
44 #include <tina/math/mathPro.h>
45 #include <tina/geometry/geomDef.h>
46 #include <tina/geometry/geomPro.h>
47 #include <tina/image/imgDef.h>
48 #include <tina/image/imgPro.h>
49 #include <tina/vision/visDef.h>
50 #include <tina/vision/visPro.h>
51
52 #include <tina/vision/vis_CorrDef.h>
53 #include <tina/vision/vis_CorrPro.h>
54
55 /*
56 #define VERT atan2(1.0, 0.0)
57 */
58 #define ANGERR 0.85
59 static Imregion im_roi;
60 static int count, pi, pj;
61
62 void setup_corLUTs(float **alpha, float **beta, int **kappa, int compress,
63 int stretch, int width)
64 {
65 int i, j;
66 float sub_pix, step;
67
68 for (i=-compress; i<stretch+1; i++)
69 {
70 sub_pix = -((float)(width/2)) - (float)i;
71 step = 1.0 + 2.0 * (float)i / ((float)width);
72 for (j=-width/2; j<width/2; j++)
73 {
74 kappa[i][j] = tina_int(sub_pix);
75 alpha[i][j] = 1.0 - (sub_pix - (float)kappa[i][j]);
76 beta [i][j] = 1.0 - alpha[i][j];
77 sub_pix += step;
78 }
79 }
80 }
81
82 static float IM_CTF_FLOAT(Imrect *im, int y, int x, int ctf_scale)
83 {
84 int i, j;
85 float sum = 0;
86
87 y *= ctf_scale;
88 x *= ctf_scale;
89
90 for (i = y; i < y + ctf_scale; i++)
91 for (j = x; j < x + ctf_scale; j++)
92 sum += IM_FLOAT(im, i, j);
93
94 return(sum / (float)(ctf_scale * ctf_scale));
95 }
96
97 float **auto_correlate(Imrect *im, Imregion *roi, int width, int compress,
98 int stretch, float *auto_partial, float *adj_partial,
99 float **norm, int **kappa, float **alpha, float **beta,
100 int ctf_scale)
101 {
102 int lx, ux, ly, uy;
103 int i, j, k, l, width_over2, j_plus1;
104 int stretch_plus1, upper_lim;
105 float sum, adj_sum, pix, a, b;
106 float *col, *adj_col, *swap_temp;
107 int *kappa_row;
108 float *alpha_row, *beta_row;
109
110 lx = roi->lx;
111 ly = roi->ly;
112 ux = roi->ux;
113 uy = roi->uy;
114 width_over2 = width/2;
115 upper_lim = ux-stretch-width_over2 + 1;
116 stretch_plus1 = stretch+1;
117
118 col = fvector_alloc(ly, uy);
119 adj_col = fvector_alloc(ly, uy);
120
121 for (k=ly; k<uy; k++)
122 {
123 adj_col[k] = IM_CTF_FLOAT(im, k, lx, ctf_scale);
124 }
125 for (j=lx; j<ux; j++)
126 {
127 swap_temp = col;
128 col = adj_col;
129 adj_col = swap_temp;
130 j_plus1 = j+1;
131
132 sum = adj_sum = 0;
133 for (k=ly; k<uy; k++)
134 {
135 pix = col[k];
136 sum += pix * pix;
137 adj_sum += pix * (adj_col[k] = IM_CTF_FLOAT(im, k, j_plus1, ctf_scale));
138 }
139 auto_partial[j] = sum;
140 adj_partial[j] = adj_sum;
141 }
142
143 for (l=lx+stretch+width_over2; l<upper_lim; l++)
144 {
145 for (i=-compress; i<stretch_plus1; i++)
146 {
147 kappa_row = kappa[i];
148 alpha_row = alpha[i];
149 beta_row = beta[i];
150 sum = 0;
151 for (j=-width_over2; j<width_over2; j++)
152 {
153 k = kappa_row[j];
154 a = alpha_row[j];
155 b = beta_row[j];
156
157 sum += a*a*auto_partial[l+k] + 2.0*a*b*adj_partial[l+k] +
158 b*b*auto_partial[l+k+1];
159 }
160 norm[l][i] = sum;
161 }
162 }
163 return(norm);
164 }
165
166 float **partial_sum(Imrect *im1, Imrect *im2, Imregion *roi1, Imregion *roi2,
167 Imregion *disp, int width, int stretch, int compress,
168 float **cross_partial, float *norm1, int ctf_scale)
169 {
170 int lx1, ux1, ly, uy;
171 int lx2, ux2;
172 int i, j, l, d, row_offset;
173 int disp_min;
174 float sum, auto1_sum, pix;
175 float *col;
176
177 lx1 = roi1->lx;
178 ly = roi1->ly;
179 ux1 = roi1->ux;
180 uy = roi1->uy;
181 lx2 = roi2->lx;
182 ux2 = roi2->ux;
183
184 disp_min = disp->lx;
185 row_offset = lx1 + width/2;
186 col = fvector_alloc(ly, uy);
187
188 auto1_sum = 0;
189 for (j=lx1; j<ux1; j++)
190 {
191 for (i=ly; i<uy; i++)
192 {
193 pix = col[i] = IM_CTF_FLOAT(im1, i, j, ctf_scale);
194 auto1_sum += pix*pix;
195 }
196
197 for (l=lx2, d=disp_min-width/2-stretch; l<ux2; l++, d++)
198 {
199 sum = 0;
200 for (i=ly; i<uy; i++)
201 {
202 sum += col[i] * IM_CTF_FLOAT(im2, i, l, ctf_scale);
203 }
204 cross_partial[j-row_offset][d] = sum;
205 }
206 }
207 fvector_free(col, ly);
208
209 *norm1 = auto1_sum;
210 return(cross_partial);
211 }
212
213 float **partial_accumulate(float **cross_partial, int **kappa, float **alpha,
214 float **beta, int stretch, int compress, int width,
215 Imregion *disp_lim, float **scores,
216 float norm1, float **norm2)
217 {
218 int i, j, d, disp_max, disp_min;
219 int width_over2, stretch_plus1;
220 int *kappa_row, k;
221 float *alpha_row, *beta_row, *scores_row, *norm2_row, sum;
222
223 disp_min = disp_lim->lx;
224 disp_max = disp_lim->ux;
225 stretch_plus1 = stretch+1;
226 width_over2 = width/2;
227
228 for (d=disp_min; d<=disp_max; d++)
229 {
230 for (i=-compress; i<stretch_plus1; i++)
231 {
232 kappa_row = kappa[i];
233 alpha_row = alpha[i];
234 beta_row = beta[i];
235 sum = 0;
236 for (j=-width_over2; j<width_over2; j++)
237 {
238 k = kappa_row[j];
239 sum += alpha_row[j] * cross_partial[j][d+k] +
240 beta_row[j] * cross_partial[j][d+k+1];
241 }
242 scores[d][i] = sum;
243 }
244 }
245
246 for (d=disp_min; d<=disp_max; d++)
247 {
248 scores_row = scores[d];
249 norm2_row = norm2[d];
250 for (j=-compress; j<stretch_plus1; j++)
251 {
252 scores_row[j] /= sqrt(norm1 * norm2_row[j]);
253 }
254 }
255
256 return(scores);
257 }
258
259 int im_count_edges(Imrect *im, Imregion *roi)
260 {
261 int lx, ly, ux, uy, i, j;
262 Imregion *valid_region=NULL;
263 int count=0;
264
265 valid_region = roi_inter(roi,im->region);
266
267 if (valid_region == NULL)
268 {
269 error("No valid region found \n", non_fatal);
270 return(0);
271 }
272 lx = valid_region->lx;
273 ly = valid_region->ly;
274 ux = valid_region->ux;
275 uy = valid_region->uy;
276
277 for (i = ly; i < uy; i++)
278 for (j = lx; j < ux; j++)
279 {
280 if (IM_PTR(im, i, j))
281 count++;
282 }
283
284 rfree(valid_region);
285 return(count);
286 }
287
288
289 static void update_maxmin(Imrect *im, float *max, float *min)
290 {
291 float pix;
292
293 if (pi>=im_roi.ly && pi<im_roi.uy && pj>=im_roi.lx && pj<im_roi.ux)
294 {
295 pix = IM_FLOAT(im, pi, pj);
296 if (pix != NO_DISP)
297 {
298 if (pix > *max) *max = pix;
299 if (pix < *min) *min = pix;
300 count++;
301 }
302 }
303 }
304
305 int imf_nz_maxmin(Imrect *im, float *max, float *min, int cen_i, int cen_j,
306 int srch_range)
307 {
308 int range;
309 float pix;
310
311 *max = -FLT_MAX;
312 *min = FLT_MAX;
313 count = 0;
314
315 im_roi.lx = im->region->lx;
316 im_roi.ux = im->region->ux;
317 im_roi.ly = im->region->ly;
318 im_roi.uy = im->region->uy;
319
320 pix = IM_FLOAT(im, cen_i, cen_j);
321 if (pix != NO_DISP)
322 {
323 if (pix > *max) *max = pix;
324 if (pix < *min) *min = pix;
325 count++;
326 }
327
328 for (range = 1; range <= srch_range; range++)
329 {
330 for (pi=cen_i-range, pj=cen_j-range; pi<=cen_i+range; pi++)
331 update_maxmin(im, max, min);
332 for (pi=cen_i-range, pj=cen_j+range; pi<=cen_i+range; pi++)
333 update_maxmin(im, max, min);
334 for (pi=cen_i-range, pj=cen_j-range+1; pj<cen_j+range; pj++)
335 update_maxmin(im, max, min);
336 for (pi=cen_i+range, pj=cen_j-range+1; pj<cen_j+range; pj++)
337 update_maxmin(im, max, min);
338 }
339
340 return(count);
341 }
342
343 static double edgel_realign(Edgel *algn_pix, Edgel *ref_pix)
344 {
345 double deltay;
346 double pos;
347
348 if (ref_pix == NULL)
349 return(vec2_x(algn_pix->pos));
350
351 deltay = vec2_y(ref_pix->pos) - vec2_y(algn_pix->pos);
352 pos = deltay/tan(algn_pix->orient);
353
354 return(vec2_x(algn_pix->pos)-pos);
355 }
356
357 float imf_hnonzero(Imrect *im, int i, int j, int range)
358 {
359 int lx, ux, ly, uy, l, col;
360 double pix;
361
362 lx = im->region->lx;
363 ux = im->region->ux;
364 ly = im->region->ly;
365 uy = im->region->uy;
366 if (i < ly || i >= uy || j < lx || j >= ux)
367 return(NO_DISP);
368
369 pix = IM_FLOAT(im, i, j);
370 if (pix != NO_DISP)
371 return(pix);
372
373 for (l=1; l<= range; l++)
374 {
375 col = j-l;
376 if (col >= lx)
377 {
378 pix = IM_FLOAT(im, i, col);
379 if (pix != NO_DISP)
380 return(pix);
381 }
382
383 col = j+l;
384 if (col < ux)
385 {
386 pix = IM_FLOAT(im, i, col);
387 if (pix != NO_DISP)
388 return(pix);
389 }
390 }
391 return(NO_DISP);
392 }
393
394 static Edgel *edge_im_hnonzero(Imrect *edge_im, int y, int x, double maxdist)
395 {
396 Edgel *edge = NULL, *best = NULL;
397 double disp, best_pos = FLT_MAX;
398 int i, lx, ux, ly, uy;
399
400 if (edge_im->vtype != ptr_v)
401 return(NULL);
402
403 lx = edge_im->region->lx;
404 ux = edge_im->region->ux;
405 ly = edge_im->region->ly;
406 uy = edge_im->region->uy;
407 if (x < lx || x >= ux || y < ly || y >= uy)
408 return(NULL);
409
410 for (i = x - tina_rint(maxdist); i <= x + tina_rint(maxdist); i++)
411 {
412 if (i < lx || i >= ux || y < ly || y >= uy)
413 continue;
414 if ((edge = IM_PTR(edge_im, y, i)) == NULL)
415 continue;
416
417 disp = fabs(vec2_x(edge->pos)-(double)x);
418 if ((disp < best_pos) && (disp < maxdist || !maxdist))
419 {
420 best = edge;
421 best_pos = disp;
422 }
423 }
424
425 return(best);
426 }
427
428 void stretch_output_block(Imrect *out_disp_im, Imrect *im1_edges,
429 Imrect *im2_edges, Imregion *roi1,
430 Imregion *roi2, double edge_delta)
431 {
432 int i, j1, j2, k, l;
433 int height, width1, width2;
434 float roi2_pos;
435 float estimate;
436 double pos1, pos2;
437 double orr_diff, xscale;
438 Edgel *edge1, *edge2;
439 Match *match;
440 /*
441 shistogram **hists;
442 hists = (shistogram **)hist_vec();
443 */
444
445
446 height = roi1->uy - roi1->ly;
447 width1 = roi1->ux - roi1->lx;
448 width2 = roi2->ux - roi2->lx;
449 i = roi1->ly + (height / 2);
450 j1 = roi1->lx + (width1 / 2);
451 j2 = roi2->lx + (width2 / 2);
452 xscale = (float)width2 / (float)width1;
453
454 for (k = -height/2; k < height/2; k++)
455 {
456 for (l = -width1/2; l < width1/2; l++)
457 {
458 /* calculate the edge disparity */
459 if ((edge1 = edge_im_hnonzero(im1_edges, i+k, j1+l, 0.0)))
460 {
461 pos1 = vec2_x(edge1->pos) - (double)j1;
462 roi2_pos = (float)j2 + ((float)pos1 * xscale);
463
464 if ((edge2 = edge_im_hnonzero(im2_edges, i+k, (int)roi2_pos, edge_delta)))
465 {
466 orr_diff = edge1->orient - atan(xscale*tan(edge2->orient));
467 if (orr_diff>PI) orr_diff -= PI;
468 if (orr_diff<-PI) orr_diff += PI;
469 if (orr_diff> PI/2.0) orr_diff = orr_diff - PI;
470 if (orr_diff< -PI/2.0) orr_diff = orr_diff + PI;
471
472 /*
473 hfill1(hists[0],orr_diff,1.0);
474 */
475
476 pos2 = edgel_realign(edge2, edge1);
477 estimate = pos2 - ((double)j1 + pos1);
478
479 /* add match property to allow uniqueness checking
480 check for existing match from previous run first */
481 if (prop_present(edge1->props, MATCH))
482 edge1->props = proplist_rm(edge1->props, MATCH);
483
484 if (fabs(orr_diff)< ANGERR)
485 {
486 IM_FLOAT(out_disp_im, i+k, j1+l) = estimate;
487 match = match_make((void *)edge1, (void *)edge2, MATCH);
488 edge1->props = proplist_addifnp(edge1->props, match, MATCH,
489 match_free, true);
490 }
491 }
492 }
493 }
494 }
495 }
496
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.