~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/image/imgColour_segment.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.