1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File :
37 * Date :
38 * Version :
39 * CVS Id :
40 *
41 * Notes :
42 *
43 * Colour segmentation functions: see Tina memo 2001-015 for an explanation.
44 *
45 *********
46 */
47
48 #include "imgColour_segment.h"
49
50 #if HAVE_CONFIG_H
51 # include <config.h>
52 #endif
53
54 #include<math.h>
55 #include<float.h>
56 #include<tina/math/mathPro.h>
57 #include<tina/image/img_GenDef.h>
58 #include<tina/image/img_GenPro.h>
59
60
61 #include "imgColour_convert.h"
62
63 #undef MAXFLOAT /* in math.h on Solaris */
64 #define MAXFLOAT FLT_MAX
65 #define MINFLOAT FLT_MIN
66 #define KNNTYPE 1000
67 #define NEIGHBOURS 20
68
69 static int no_peaks=0;
70 static int no_nodes=0;
71 static int abort_flag=0;
72 static int total_knn=0;
73 static double pi = 3.1415926535;
74 List *density_list = NULL;
75
76
77 static void dot_product(float *dot_p, float *data1, float *data2, float *data3, int vectordim)
78 {
79 int i;
80
81 *dot_p = 0;
82 for(i=0; i<vectordim; i++)
83 {
84 *dot_p+= ((data2[i]/data2[(vectordim+i)])-(data1[i]/data1[(vectordim+i)]))
85 *((data3[i]/data3[(vectordim+i)])-(data1[i]/data1[(vectordim+i)]));
86 }
87 }
88
89
90 static void magnitude(float *mag, float *data1, float *data3, int vectordim)
91 {
92 int i;
93
94 *mag = 0;
95 for(i=0; i<vectordim; i++)
96 {
97 *mag += sqr((data3[i] / data3[(vectordim+i)]) - (data1[i] / data1[(vectordim+i)]));
98 }
99 *mag = sqrt(*mag);
100 }
101
102
103 static double distance_func(float *data1, float *data2, int n)
104 {
105 int i;
106 double dist = 0.0;
107
108 for (i=0;i<n;i++)
109 {
110 dist += sqr( ((data1[i] / data1[(n+i)]) - (data2[i] / data2[(n+i)])) );
111 }
112 return(dist/2.0);
113 }
114
115
116 static double sortfunc(Match *element, int type)
117 {
118 double val=0.0;
119
120 val = (double)(element->weight);
121 return (-val);
122 }
123
124
125 static void **knn_listfree(void **nearest, int near)
126 {
127 int i;
128
129 for (i=0;i<near;i++) rfree((void *)(nearest[i]));
130 pvector_free(nearest,0);
131 return(NULL);
132 }
133
134
135 static void **knn_list(List *density_list, int neighbours, float *data,
136 double(*distance_func)(/* */))
137 {
138 List *plist=NULL;
139 void **nearest=NULL;
140 Match *mptr=NULL, *nptr=NULL, *lptr=NULL;
141 float knn_dist=MAXFLOAT;
142 double dist;
143 int i;
144
145 if (density_list==NULL) return(NULL);
146 nearest = pvector_alloc(0,neighbours);
147 for(plist = density_list; plist !=NULL; plist = plist->next)
148 {
149 mptr = plist->to;
150 dist = distance_func(data, mptr->to1, mptr->type);
151
152 /* now convert to a Euclidean distance metric */
153 dist = sqrt((2*dist));
154 if (dist < knn_dist)
155 {
156 nptr = NULL;
157 nptr = match_alloc(mptr->type);
158 nptr->to1 = data;
159 nptr->weight = dist;
160 nptr->to2 = mptr; /* store matching node on this match */
161 for (i=0; i<neighbours; i++)
162 {
163 lptr = (Match *)nearest[i];
164 if (nptr!=NULL && (lptr ==NULL || dist < lptr->weight))
165 {
166 nearest[i] = nptr;
167 if ((nptr!=NULL)&&(i==neighbours-1))
168 knn_dist = nptr->weight;
169 nptr = lptr;
170 }
171 }
172 if (nptr!=NULL) rfree(nptr);
173 }
174 }
175 return(nearest);
176 }
177
178
179 List *knn_destroy(List *density_list)
180 {
181 List *nlist=NULL;
182 Match *nptr=NULL;
183 float *data=NULL;
184
185 for (nlist = density_list; nlist!=NULL; nlist = nlist->next)
186 {
187 nptr = nlist->to;
188 data = nptr->to1;
189 fvector_free(data,0);
190 match_free(nptr);
191 rfree((void *)nlist);
192 }
193 total_knn = 0;
194 density_list=NULL;
195 return(NULL);
196 }
197
198
199 static List *peak_list_destroy(List *peak_list)
200 {
201 List *nlist=NULL;
202 Match *nptr=NULL;
203
204 for (nlist = peak_list; nlist!=NULL; nlist = nlist->next)
205 {
206 nptr = nlist->to;
207 match_free((void *)nptr);
208 rfree((void *)nlist);
209 }
210 return(NULL);
211 }
212
213
214 static List *knn_build(List *density_list, float *data, void *classdata,
215 double(*distance_func)(/* */), float threshold, int vectordim)
216 {
217 void **nearest = NULL;
218 Match *nptr=NULL;
219 float *data_copy=NULL;
220 int i;
221
222 nearest = knn_list(density_list, 1, data, distance_func);
223 if(nearest!=NULL && nearest[0]==NULL)
224 {
225 abort_flag=1;
226 error("Internal error type 1: see documentation", warning);
227 return(density_list);
228 }
229 if (nearest==NULL || ((Match *)(nearest[0]))->weight > threshold)
230 {
231 nptr = match_alloc(vectordim);
232 data_copy = (float *)fvector_alloc(0, (2*vectordim));
233 for (i=0;i<(2*vectordim);i++) data_copy[i] = data[i];
234 nptr->to1 = data_copy;
235 nptr->to2 = classdata;
236 density_list = link_addtostart(density_list, link_alloc(nptr, KNNTYPE));
237 no_nodes = no_nodes+1;
238 total_knn++;
239 }
240 if (nearest!=NULL) nearest = knn_listfree(nearest, 1);
241 return(density_list);
242 }
243
244
245 static void knn_fill(List *density_list, float *data, double(*distance_func)(/* */), int vectordim)
246 {
247 int i;
248 void **nearest=NULL;
249 Match *mptr=NULL, *nptr=NULL;
250 float dist, extra;
251 int neighbours = NEIGHBOURS;
252
253 nearest = knn_list(density_list, neighbours, data, distance_func);
254 for(i=0; i<neighbours; i++)
255 {
256 if((mptr = nearest[i]) == NULL) break;
257 nptr = mptr->to2;
258 dist = mptr->weight;
259 extra = exp(-dist);
260 nptr->weight += extra; /* accumulate density of patterns */
261 }
262 nearest = knn_listfree(nearest,neighbours);
263 }
264
265
266 static List *knn_fill_2v(List *density_list, double (*distance_func)(/* */),
267 int vectordim, double threshold)
268 {
269 List *nlist=NULL, *mlist=NULL;
270 Match *mptr=NULL, *nptr=NULL, *lptr=NULL, *optr=NULL, *topptr=NULL;
271 float *data1=NULL, *data2=NULL, *data3=NULL;
272 void **nearest=NULL;
273 int i, neighbours, voronoid_flag=0;
274 double dist, gradient, best_gradient=MINFLOAT, height_there, height_here;
275 float voronoid_dist, dot_p, mag;
276
277 neighbours = no_nodes;
278 for(nlist = density_list; nlist!=NULL; nlist = nlist->next)
279 {
280 topptr = NULL;
281 best_gradient = -MAXFLOAT;
282 mptr = nlist->to;
283 data1 = mptr->to1;
284 height_here = mptr->weight;
285 nearest = knn_list(density_list, no_nodes, data1, distance_func);
286 for (i=0; i<neighbours; i++)
287 {
288 nptr = (Match *)nearest[i];
289 if(nptr==NULL) break;
290 dist = sqrt(nptr->weight);
291 lptr = nptr->to2;
292 data2 = lptr->to1;
293
294 height_there = lptr->weight;
295 gradient = (height_there/(sqrt(2*pi)))
296 * exp(-pow(dist, 2)/2);
297 voronoid_flag=0;
298
299 if(dist!=0&&nptr!=NULL&&gradient>best_gradient&&height_here<height_there)
300 {
301 for(mlist = density_list; mlist!=NULL; mlist=mlist->next)
302 {
303 optr=mlist->to;
304 data3 = optr->to1;
305 dot_product(&dot_p, data1, data2, data3, vectordim);
306 magnitude(&mag, data1, data3, vectordim);
307 voronoid_dist = sqr(mag) / dot_p;
308
309 if(voronoid_dist<1.0&&voronoid_dist>0&&data3!=data2)
310 {
311 voronoid_flag=1;
312 }
313 }
314 if(voronoid_flag==0)
315 {
316 best_gradient = gradient;
317 topptr = (Match *)lptr;
318 }
319 }
320 }
321 mptr->to2 = topptr;
322 }
323 nearest = knn_listfree(nearest, no_nodes);
324 return density_list;
325 }
326
327
328 static List *hill_climb(List *density_list, float *data, double(*distance_func)(/* */),
329 List *peak_list, int vectordim, double threshold)
330 {
331 /* This subroutine performs the hill climb to return a label for the data point */
332
333 Match *mptr=NULL, *nptr=NULL, *optr=NULL, *pptr=NULL, *qptr=NULL, *rptr=NULL;
334 void **nearest=NULL, *topptr=NULL;
335 List *nlist=NULL, *plist=NULL;
336 float *data2=NULL, *data3=NULL;
337 double dist, best_gradient = -MINFLOAT, gradient, height_there;
338 int neighbours, i, voronoid_flag=0;
339 float voronoid_dist, dot_p, mag;
340
341 /* First find the nearest neighbour to the current data point from density list */
342
343 neighbours = no_nodes;
344
345 nearest = knn_list(density_list, neighbours, data, distance_func);
346 for(i=0; i<neighbours; i++)
347 {
348 mptr = (Match *)nearest[i];
349 if(mptr==NULL) break;
350 dist = sqrt(mptr->weight);
351 nptr = (Match *)mptr->to2;
352 data2 = nptr->to1;
353 height_there = nptr->weight;
354 gradient = (height_there/(sqrt(2*pi))) * exp(-pow(dist, 2)/2);
355 voronoid_flag=0;
356
357 if(dist!=0&&mptr!=NULL&&gradient>best_gradient)
358 {
359 for(plist=density_list; plist!=NULL; plist=plist->next)
360 {
361 rptr=plist->to;
362 data3 = rptr->to1;
363 dot_product(&dot_p, data, data2, data3, vectordim);
364 magnitude(&mag, data, data3, vectordim);
365 voronoid_dist = sqr(mag) / dot_p;
366 if(voronoid_dist<1.0&&voronoid_dist>0.0&&data3!=data2)
367 {
368 voronoid_flag=1;
369 }
370 }
371 if(voronoid_flag==0)
372 {
373 best_gradient = gradient;
374 topptr = (Match *)nptr;
375 }
376 }
377 }
378 optr = topptr;
379
380 /* Now climb the hill */
381
382 while(optr!=NULL)
383 {
384 nptr = optr;
385 optr = nptr->to2;
386 }
387
388 /* At this point nptr should point to the element of density list that lies at
389 the top of the hill. Now compile the peak list and accumulate the contributions
390 to each peak. */
391
392 for(nlist=peak_list; nlist!=NULL; nlist=nlist->next)
393 {
394 qptr = ((Match *)(nlist->to));
395 if(nptr==(qptr->to1))
396 {
397 (qptr->weight)=(qptr->weight)+1;
398 nearest = knn_listfree(nearest, neighbours);
399 return peak_list;
400 }
401 }
402
403 pptr = match_alloc(vectordim);
404 pptr->to1 = nptr;
405 pptr->to2 = NULL;
406 pptr->weight = 1;
407 peak_list = link_addtostart(peak_list, link_alloc(pptr, KNNTYPE));
408 no_peaks = no_peaks+1;
409 nearest = knn_listfree(nearest, neighbours);
410 return peak_list;
411 }
412
413
414 static List *peak_list_sort(List *peak_list, int vectordim)
415 {
416
417 /* Function to sort the peak list and create the histogram */
418
419 List *peak_list_2=NULL, *nlist=NULL;
420 Match *lptr=NULL;
421 int i=0;
422
423 for(nlist=peak_list; nlist!=NULL; nlist=nlist->next)
424 {
425 i=i+1;
426 lptr = (Match *)nlist->to;
427 peak_list_2 = list_add_sorted(peak_list_2,
428 link_alloc(lptr, KNNTYPE), sortfunc);
429 }
430
431 list_rm_links(peak_list);
432 peak_list = NULL;
433 return peak_list_2;
434 }
435
436
437 static int label_assign(List *density_list, float *data, double(*distance_func)(/* */),
438 List *peak_list, int vectordim)
439 {
440 int label, i;
441 Match *mptr=NULL, *nptr=NULL, *optr=NULL, *qptr=NULL;
442 void **nearest=NULL;
443 List *nlist=NULL;
444
445 nearest = knn_list(density_list, 1, data, distance_func);
446 mptr = (Match *)nearest[0];
447 nptr = mptr->to2;
448 optr = nptr->to2; /* This is the pointer to the next density list position up the hill */
449
450 while(optr!=NULL)
451 {
452 nptr = optr;
453 optr = nptr->to2;
454 }
455
456 i=0;
457 for(nlist=peak_list; nlist!=NULL; nlist=nlist->next)
458 {
459 i=i+1;
460 qptr = (Match *)(nlist->to);
461 if(nptr==(qptr->to1))
462 {
463 label = i;
464 nearest = knn_listfree(nearest, 1);
465 return label;
466 }
467 }
468
469 label = 1.0;
470 return label;
471 /* ?? not reached */
472 nearest = knn_listfree(nearest,1);
473 }
474
475
476 Match *colour_segment(Imrect *im_j, Imrect *im_k, void ***covar_matrix, double threshold)
477 {
478 Match *return_match=NULL;
479 Imrect *im = NULL;
480 Imregion *roi=NULL;
481 float *row1=NULL, *data=NULL, *covar_ptr=NULL;
482 int i, j;
483 List *peak_list = NULL;
484
485 no_peaks = 0;
486 no_nodes = 0;
487
488 if(density_list!=NULL)
489 {
490 density_list = knn_destroy(density_list);
491 }
492
493 roi = im_j->region;
494 data = (float *)fvector_alloc(0, 4);
495 im = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), roi, float_v);
496 row1 = fvector_alloc(roi->lx, roi->ux);
497
498 for (i = roi->ly; i < roi->uy; i++)
499 {
500 for (j = roi->lx; j < roi->ux; j++)
501 {
502 data[0] = ((float **)im_j->data)[i][j];
503 data[1] = ((float **)im_k->data)[i][j];
504 covar_ptr = (float *)covar_matrix[i][j];
505 data[2] = covar_ptr[0];
506 data[3] = covar_ptr[1];
507
508 density_list = knn_build(density_list,data, NULL, distance_func,
509 threshold, 2);
510 if(abort_flag==1)
511 {
512 error("Leaving segmentation", warning);
513 im_free(im);
514 return(NULL);
515 }
516 }
517 }
518 for (i = roi->ly; i < roi->uy; i++)
519 {
520 for (j = roi->lx; j < roi->ux; j++)
521 {
522 data[0] = ((float **)im_j->data)[i][j];
523 data[1] = ((float **)im_k->data)[i][j];
524 covar_ptr = (float *)covar_matrix[i][j];
525 data[2] = covar_ptr[0];
526 data[3] = covar_ptr[1];
527
528 knn_fill(density_list,data, distance_func, 2);
529 }
530 }
531
532 density_list = knn_fill_2v(density_list, distance_func,
533 2, threshold);
534
535 for (i = roi->ly; i < roi->uy; i++)
536 {
537 for (j = roi->lx; j < roi->ux; j++)
538 {
539 data[0] = ((float **)im_j->data)[i][j];
540 data[1] = ((float **)im_k->data)[i][j];
541 covar_ptr = (float *)covar_matrix[i][j];
542 data[2] = covar_ptr[0];
543 data[3] = covar_ptr[1];
544
545 peak_list = hill_climb(density_list,data,
546 distance_func, peak_list, 2, threshold);
547 }
548 }
549 peak_list = peak_list_sort(peak_list, 2);
550
551 for (i = roi->ly; i < roi->uy; i++)
552 {
553 for (j = roi->lx; j < roi->ux; j++)
554 {
555 data[0] = ((float **)im_j->data)[i][j];
556 data[1] = ((float **)im_k->data)[i][j];
557 covar_ptr = (float *)covar_matrix[i][j];
558 data[2] = covar_ptr[0];
559 data[3] = covar_ptr[1];
560
561 row1[j] = (float)label_assign(density_list,data,
562 distance_func, peak_list, 2);
563 }
564 im_put_rowf(row1, im, i, roi->lx, roi->ux);
565 }
566
567 format("The number of nodes is %d\n", no_nodes);
568 format("The number of peaks is %d\n", no_peaks);
569
570 fvector_free(row1, roi->lx);
571 fvector_free(data, 0);
572 peak_list = peak_list_destroy(peak_list);
573 return_match = match_alloc(1010);
574 return_match->to1 = im;
575 return_match->to2 = density_list;
576 density_list=NULL; /* Will be freed in the claaing function */
577 return(return_match);
578 }
579
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.