segmentation
All Data Structures Namespaces Files Functions Variables Modules Pages
BgEdgeDetect.cpp
1 // Name: BgEdgeDetect.cpp
3 // Purpose: BgEdgeDetect class functions
4 // Author: Bogdan Georgescu
5 // Modified by:
6 // Created: 06/22/2000
7 // Copyright: (c) Bogdan Georgescu
8 // Version: v0.1
10 
11 #include <math.h>
12 #include "edge/BgImage.h"
13 #include "edge/BgEdge.h"
14 #include "edge/BgEdgeList.h"
15 #include "edge/BgEdgeDetect.h"
16 #include "edge/BgDefaults.h"
17 #include <stdio.h>
18 
19 #define TOL_E 2.2e-8
20 //#define TOL_E 0.05
21 #define SIGN(x) (x<0)? -1:1;
22 
23 static int my_sign(double val)
24 {
25  if(val>TOL_E)
26  return 1;
27  if(val<-TOL_E)
28  return -1;
29  return 0;
30 }
31 
32 extern double factorial(double);
33 
34 BgEdgeDetect::BgEdgeDetect(int filtDim)
35 {
36  havePerm_ = false;
37  WL_ = filtDim;
38  WW_ = 2*WL_+1;
39  nhcust_ = 0;
40  nlcust_ = 0;
41  tcustx_ = new float[MAX_CUSTT];
42  tcusty_ = new float[MAX_CUSTT];
43  CreateFilters();
44  CreateLookTable();
45 }
46 
47 BgEdgeDetect::~BgEdgeDetect()
48 {
49  if (havePerm_==true)
50  {
51  delete [] permGx_;
52  delete [] permGy_;
53  delete [] permConf_;
54  delete [] permRank_;
55  delete [] permNmxRank_;
56  delete [] permNmxConf_;
57  }
58  if (nhcust_>0)
59  {
60  delete [] hcustx_;
61  delete [] hcusty_;
62  }
63  if (nlcust_>0)
64  {
65  delete [] lcustx_;
66  delete [] lcusty_;
67  }
68  delete [] tcustx_;
69  delete [] tcusty_;
70  DeleteLookTable();
71 }
72 
73 void BgEdgeDetect::IsGood(void)
74 {
75  if (havePerm_==true)
76  bgLog("good\n");
77  else
78  bgLog("bad\n");
79 }
80 
81 
82 float BgEdgeDetect::EllipseEval(float x, float y)
83 {
84  return ((x*x)/(rankTr_*rankTr_)+(y*y)/(confTr_*confTr_)-1);
85 }
86 
87 float BgEdgeDetect::EllipseComp(float x0, float y0, float x, float y)
88 {
89 // return (EllipseEval(x,y)-EllipseEval(x0,y0));
90  return ((x*x-x0*x0)/(rankTr_*rankTr_)+(y*y-y0*y0)/(confTr_*confTr_));
91 }
92 
93 float BgEdgeDetect::LineEval(float x, float y)
94 {
95  return (confTr_*x+rankTr_*y-confTr_*rankTr_);
96 }
97 
98 float BgEdgeDetect::LineComp(float x0, float y0, float x, float y)
99 {
100 // return (LineEval(x,y)-LineEval(x0,y0));
101  return (confTr_*(x-x0)+rankTr_*(y-y0));
102 }
103 
104 float BgEdgeDetect::VerticalLineEval(float x, float y)
105 {
106  return (x-rankTr_);
107 }
108 
109 float BgEdgeDetect::VerticalLineComp(float x0, float y0, float x, float y)
110 {
111 // return (VerticalLineEval(x,y)-VerticalLineEval(x0,y0));
112  return (x-x0);
113 
114 }
115 
116 float BgEdgeDetect::HorizontalLineEval(float x, float y)
117 {
118  return(y-confTr_);
119 }
120 
121 float BgEdgeDetect::HorizontalLineComp(float x0, float y0, float x, float y)
122 {
123 // return (HorizontalLineEval(x,y)-HorizontalLineEval(x0,y0));
124  return (y-y0);
125 }
126 
127 float BgEdgeDetect::SquareEval(float x, float y)
128 {
129  if ((x/rankTr_)>(y/confTr_))
130  return(x-rankTr_);
131  else
132  return(y-confTr_);
133 }
134 
135 float BgEdgeDetect::SquareComp(float x0, float y0, float x, float y)
136 {
137 // return(SquareEval(x,y)-SquareEval(x0,y0));
138  static float tret;
139  tret = ((x/rankTr_)>(y/confTr_)) ? x-rankTr_ : y-confTr_;
140  tret -= ((x0/rankTr_)>(y0/confTr_)) ? x0-rankTr_ : y0-confTr_;
141  return tret;
142 }
143 
144 float BgEdgeDetect::CustomRegionEval(float r,float c)
145 {
146  //evaluate user region function
147  //returns -1 if inside +1 if outside
148 
149  if ((r+c)<=ZERO_TRESH)
150  return -1;
151  int i;
152  int crossings=0;
153  float x;
154 
155  //shift to origin
156  for (i=0; i<ncust_; i++)
157  {
158  tcustx_[i]=custx_[i]-r;
159  tcusty_[i]=custy_[i]-c;
160  }
161 
162  for (i=0; i<(ncust_-1); i++)
163  {
164  if ( (tcusty_[i] >0 && tcusty_[i+1]<=0) ||
165  (tcusty_[i+1]>0 && tcusty_[i] <=0) )
166  {
167  x = (tcustx_[i]*tcusty_[i+1]-tcustx_[i+1]*tcusty_[i])/(tcusty_[i+1]-tcusty_[i]);
168  if (x>0)
169  crossings++;
170  }
171  }
172 
173  if ((crossings % 2) ==1)
174  return -1;
175  else
176  return 1;
177 }
178 
179 float BgEdgeDetect::CustomRegionComp(float r0, float c0, float r, float c)
180 {
181  return 0;
182 }
183 
184 void BgEdgeDetect::GenerateMaskAngle(double* a,double theta) {
185  static int sflag;
186  static int i,j,k;
187  static double cval[4];
188  static double corner[2][4];
189  static double sinv,cosv;
190  static double intrs[2][4];
191  static int scor[4],nscor[4];
192  static int sind,rowind,colind;
193  static double cordi[2][4];
194  static int lsigind,corin;
195  static int sigind[4];
196  static double diffin[2];
197  static double comcoor;
198 
199  theta = theta*PI/180.0;
200  sinv = sin(theta);
201  cosv = cos(theta);
202 
203  for (i=0; i<WW_*WW_; i++)
204  a[i]=0;
205 
206  for (i=WL_; i>=-WL_; i--)
207  {
208  for(j=-WL_; j<=WL_; j++)
209  {
210  corner[0][0] = j-0.5;
211  corner[0][1] = j+0.5;
212  corner[0][2] = j+0.5;
213  corner[0][3] = j-0.5;
214 
215  corner[1][0] = i+0.5;
216  corner[1][1] = i+0.5;
217  corner[1][2] = i-0.5;
218  corner[1][3] = i-0.5;
219 
220  cval[0] = -sinv*corner[0][0]+cosv*corner[1][0];
221  cval[1] = -sinv*corner[0][1]+cosv*corner[1][1];
222  cval[2] = -sinv*corner[0][2]+cosv*corner[1][2];
223  cval[3] = -sinv*corner[0][3]+cosv*corner[1][3];
224 
225  scor[0] = my_sign(cval[0]);
226  scor[1] = my_sign(cval[1]);
227  scor[2] = my_sign(cval[2]);
228  scor[3] = my_sign(cval[3]);
229 
230  sind = 0;
231  if (scor[0]!=0)
232  nscor[sind++] = scor[0];
233  if (scor[1]!=0)
234  nscor[sind++] = scor[1];
235  if (scor[2]!=0)
236  nscor[sind++] = scor[2];
237  if (scor[3]!=0)
238  nscor[sind++] = scor[3];
239 
240  sflag = 0;
241  for (k=1;k<sind;k++)
242  {
243  if (nscor[k]!=nscor[0])
244  sflag++;
245  }
246 
247  rowind = i+WL_;
248  colind = j+WL_;
249 
250  if (sflag==0)
251  {
252  if (nscor[0]==1)
253  a[colind+rowind*WW_] = 1.0;
254  else
255  a[colind+rowind*WW_] = 0.0;
256  }
257 
258  if (sflag!=0)
259  {
260  for (k=0; k<4; k++)
261  intrs[0][k] = intrs[1][k] = 0.0;
262 
263  if(scor[0]==0)
264  {
265  intrs[0][0] = corner[0][0];
266  intrs[1][0] = corner[1][0];
267  }
268  if (scor[0]*scor[1]<0)
269  {
270  intrs[0][0] = corner[1][0]*cosv/sinv;
271  intrs[1][0] = corner[1][0];
272  }
273  if (scor[1]==0)
274  {
275  intrs[0][1] = corner[0][1];
276  intrs[1][1] = corner[1][1];
277  }
278  if (scor[1]*scor[2]<0)
279  {
280  intrs[0][1] = corner[0][1];
281  intrs[1][1] = corner[0][1]*sinv/cosv;
282  }
283  if (scor[2]==0)
284  {
285  intrs[0][2] = corner[0][2];
286  intrs[1][2] = corner[1][2];
287  }
288  if (scor[2]*scor[3]<0)
289  {
290  intrs[0][2] = corner[1][2]*cosv/sinv;
291  intrs[1][2] = corner[1][2];
292  }
293  if (scor[3]==0)
294  {
295  intrs[0][3] = corner[0][3];
296  intrs[1][3] = corner[1][3];
297  }
298  if (scor[3]*scor[0]<0)
299  {
300  intrs[0][3] = corner[0][3];
301  intrs[1][3] = corner[0][3]*sinv/cosv;
302  }
303 
304  corin = 0;
305  if (fabs(intrs[0][0])>TOL_E || fabs(intrs[1][0])>TOL_E)
306  {
307  cordi[0][corin] = intrs[0][0];
308  cordi[1][corin++] = intrs[1][0];
309  }
310  if (fabs(intrs[0][1])>TOL_E || fabs(intrs[1][1])>TOL_E)
311  {
312  cordi[0][corin] = intrs[0][1];
313  cordi[1][corin++] = intrs[1][1];
314  }
315  if (fabs(intrs[0][2])>TOL_E || fabs(intrs[1][2])>TOL_E)
316  {
317  cordi[0][corin] = intrs[0][2];
318  cordi[1][corin++] = intrs[1][2];
319  }
320  if (fabs(intrs[0][3])>TOL_E || fabs(intrs[1][3])>TOL_E)
321  {
322  cordi[0][corin] = intrs[0][3];
323  cordi[1][corin++] = intrs[1][3];
324  }
325 
326  lsigind=0;
327  if (scor[0]>0)
328  sigind[lsigind++] = 0;
329  if (scor[1]>0)
330  sigind[lsigind++] = 1;
331  if (scor[2]>0)
332  sigind[lsigind++] = 2;
333  if (scor[3]>0)
334  sigind[lsigind++] = 3;
335 
336  if (lsigind==1)
337  {
338  a[colind+rowind*WW_] = 0.5*fabs(cordi[0][0]-cordi[0][1])*fabs(cordi[1][0]-cordi[1][1]);
339  }
340  if (lsigind==2)
341  {
342  diffin[0] = (int) fabs(cordi[0][0]-cordi[0][1]);
343  diffin[1] = (int) fabs(cordi[1][0]-cordi[1][1]);
344  if (diffin[0]==1)
345  {
346  comcoor = corner[1][sigind[0]];
347  a[colind+rowind*WW_] = 0.5*(fabs(comcoor-cordi[1][0])+fabs(comcoor-cordi[1][1]));
348  }
349  if (diffin[1]==1)
350  {
351  comcoor = corner[0][sigind[0]];
352  a[colind+rowind*WW_] = 0.5*(fabs(comcoor-cordi[0][0])+fabs(comcoor-cordi[0][1]));
353  }
354  }
355  if(lsigind==3)
356  {
357  a[colind+rowind*WW_] = 1.0-0.5*fabs(cordi[0][0]-cordi[0][1])*fabs(cordi[1][0]-cordi[1][1]);
358  }
359  }
360  }
361  }
362 
363  //A=A-mean(mean(A));
364  comcoor = 0;
365  for (i=0; i<WW_*WW_; i++)
366  comcoor += a[i];
367  comcoor /= WW_*WW_;
368  for (i=0; i<WW_*WW_; i++)
369  a[i] -= comcoor;
370 
371  //A=A/norm(A,'fro')
372  comcoor = 0;
373  for (i=0; i<WW_*WW_; i++)
374  comcoor += a[i]*a[i];
375  comcoor = sqrt(comcoor);
376  for (i=0; i<WW_*WW_; i++)
377  a[i] /= comcoor;
378 }
379 
380 void BgEdgeDetect::CreateFilters(void)
381 {
382  int i,j;
383  double w;
384  for (i=-WL_; i<=WL_; i++)
385  {
386  w = pow(2.0,(-2*WL_))*factorial(2*WL_)/(factorial(WL_-i)*factorial(WL_+i));
387  smofil_[i+WL_] = w;
388  diffil_[i+WL_] = (2*i*w)/WL_;
389  }
390  for (j=0; j<WW_; j++)
391  {
392  for (i=0; i<WW_; i++)
393  {
394  wdy_[j+i*WW_] = wdx_[i+j*WW_] = smofil_[j]*diffil_[i];
395  }
396  }
397 
398  double norms = 0;
399  double normd = 0;
400  for (i=0; i<WW_; i++)
401  {
402  norms += smofil_[i]*smofil_[i];
403  normd += diffil_[i]*diffil_[i];
404  }
405 
406  for (j=0; j<WW_; j++)
407  {
408  for (i=0; i<WW_; i++)
409  {
410  mQ_[i][j] = (smofil_[j]*smofil_[i])/norms + (diffil_[j]*diffil_[i])/normd;
411  mN_[i][j] = (i==j) ? 1-mQ_[i][j] : -mQ_[i][j];
412  }
413  }
414 }
415 
416 void BgEdgeDetect::CreateLookTable()
417 {
418  bgLog("Creating angle lookup table\n");
419  int i;
420  for (i=-180; i<=180; i++)
421  {
422  lookTable_[i+180] = new double[WW_*WW_];
423  GenerateMaskAngle(lookTable_[i+180], (double) i);
424  }
425 }
426 
427 void BgEdgeDetect::DeleteLookTable()
428 {
429  int i;
430  for (i=0; i<NO_ANGLES; i++)
431  {
432  delete [] lookTable_[i];
433  }
434 }
435 
436 void BgEdgeDetect::GetPixels(int* nopix, int* pixx, int* pixy, double x1, double x2, double y1, double y2)
437 {
438  double minx,maxx,miny,maxy;
439 
440  if (x1<x2)
441  {
442  minx = x1;
443  maxx = x2;
444  }
445  else
446  {
447  minx = x2;
448  maxx = x1;
449  }
450 
451  if (y1<y2)
452  {
453  miny = y1;
454  maxy = y2;
455  }
456  else
457  {
458  miny = y2;
459  maxy = y1;
460  }
461 
462  int i,j,npix;
463  npix = 0;
464  for (j=0; j<y_; j++)
465  {
466  for (i=0; i<x_; i++)
467  {
468  if (permRank_[i+j*x_]<maxx && permRank_[i+j*x_]>minx &&
469  permConf_[i+j*x_]<maxy && permConf_[i+j*x_]>miny)
470  {
471  pixx[npix] = i;
472  pixy[npix++] = j;
473  }
474  }
475  }
476  *nopix = npix;
477 }
478 
479 void BgEdgeDetect::GetNmxPixels(int* nopix, int* pixx, int* pixy, double x1, double x2, double y1, double y2)
480 {
481  double minx,maxx,miny,maxy;
482  if (x1<x2)
483  {
484  minx = x1;
485  maxx = x2;
486  }
487  else
488  {
489  minx = x2;
490  maxx = x1;
491  }
492 
493  if (y1<y2)
494  {
495  miny = y1;
496  maxy = y2;
497  }
498  else
499  {
500  miny = y2;
501  maxy = y1;
502  }
503 
504  int i,j,npix;
505  npix = 0;
506  for (j=0; j<y_; j++)
507  {
508  for (i=0; i<x_; i++)
509  {
510  if (permNmxRank_[i+j*x_]<maxx && permNmxRank_[i+j*x_]>minx &&
511  permNmxConf_[i+j*x_]<maxy && permNmxConf_[i+j*x_]>miny)
512  {
513  pixx[npix] = i;
514  pixy[npix++] = j;
515  }
516  }
517  }
518  *nopix = npix;
519 }
520 
521 void BgEdgeDetect::DoRecompute(BgEdgeList* cel, double nmxr, double nmxc,
522  double rh, double ch, double rl, double cl,
523  int nMin, int nmxType, int hystTypeHigh, int hystTypeLow)
524 {
525  float *tr, *tc, *tdh, *tdl;
526  bgLog("Start edge detection...\n");
527  tr = new float[x_*y_];
528  tc = new float[x_*y_];
529  tdh = new float[x_*y_];
530  tdl = new float[x_*y_];
531 
532  //Nonmaximum supression
533  bgLog("...nonmaxima supression: ");
534 
535  float (BgEdgeDetect::*fcomp)(float,float,float,float);
536  float (BgEdgeDetect::*feval)(float,float);
537 
538  switch(nmxType)
539  {
540  case FC_ELLIPSE:
541  fcomp = &BgEdgeDetect::EllipseComp;
542  feval = &BgEdgeDetect::EllipseEval;
543  bgLog("arc\n");
544  break;
545  case FC_VERT_LINE:
546  fcomp = &BgEdgeDetect::VerticalLineComp;
547  feval = &BgEdgeDetect::VerticalLineEval;
548  bgLog("vertical line\n");
549  break;
550  case FC_HORIZ_LINE:
551  fcomp = &BgEdgeDetect::HorizontalLineComp;
552  feval = &BgEdgeDetect::HorizontalLineEval;
553  bgLog("horizontal line\n");
554  break;
555  case FC_SQUARE_BOX:
556  fcomp = &BgEdgeDetect::SquareComp;
557  feval = &BgEdgeDetect::SquareEval;
558  bgLog("box\n");
559  break;
560  case FC_LINE:
561  fcomp = &BgEdgeDetect::LineComp;
562  feval = &BgEdgeDetect::LineEval;
563  bgLog("line\n");
564  break;
565  default:
566  bgLog("Type not known\n");
567  return;
568  }
569  confTr_ = (float) nmxc;
570  rankTr_ = (float) nmxr;
571 
572  NewNonMaxSupress(permRank_,permConf_,permGx_,permGy_,permNmxRank_,permNmxConf_,fcomp);
573  bgLog("...hysteresis thresholding, high: ");
574 
575  switch(hystTypeHigh)
576  {
577  case FC_ELLIPSE:
578  fcomp = &BgEdgeDetect::EllipseComp;
579  feval = &BgEdgeDetect::EllipseEval;
580  bgLog("arc");
581  break;
582  case FC_VERT_LINE:
583  fcomp = &BgEdgeDetect::VerticalLineComp;
584  feval = &BgEdgeDetect::VerticalLineEval;
585  bgLog("vertical line");
586  break;
587  case FC_HORIZ_LINE:
588  fcomp = &BgEdgeDetect::HorizontalLineComp;
589  feval = &BgEdgeDetect::HorizontalLineEval;
590  bgLog("horizontal line");
591  break;
592  case FC_SQUARE_BOX:
593  fcomp = &BgEdgeDetect::SquareComp;
594  feval = &BgEdgeDetect::SquareEval;
595  bgLog("box");
596  break;
597  case FC_LINE:
598  fcomp = &BgEdgeDetect::LineComp;
599  feval = &BgEdgeDetect::LineEval;
600  bgLog("line");
601  break;
602  case FC_CUSTOM:
603  custx_ = hcustx_;
604  custy_ = hcusty_;
605  ncust_ = nhcust_;
606  fcomp = &BgEdgeDetect::CustomRegionComp;
607  feval = &BgEdgeDetect::CustomRegionEval;
608  bgLog("custom");
609  break;
610  }
611  confTr_ = (float) ch;
612  rankTr_ = (float) rh;
613  StrConfEstim(permNmxRank_, permNmxConf_, tdh, feval);
614 
615  bgLog(" low: ");
616 
617  switch(hystTypeLow) {
618  case FC_ELLIPSE:
619  fcomp = &BgEdgeDetect::EllipseComp;
620  feval = &BgEdgeDetect::EllipseEval;
621  bgLog("arc\n");
622  break;
623  case FC_VERT_LINE:
624  fcomp = &BgEdgeDetect::VerticalLineComp;
625  feval = &BgEdgeDetect::VerticalLineEval;
626  bgLog("vertical line\n");
627  break;
628  case FC_HORIZ_LINE:
629  fcomp = &BgEdgeDetect::HorizontalLineComp;
630  feval = &BgEdgeDetect::HorizontalLineEval;
631  bgLog("horizontal line\n");
632  break;
633  case FC_SQUARE_BOX:
634  fcomp = &BgEdgeDetect::SquareComp;
635  feval = &BgEdgeDetect::SquareEval;
636  bgLog("box\n");
637  break;
638  case FC_LINE:
639  fcomp = &BgEdgeDetect::LineComp;
640  feval = &BgEdgeDetect::LineEval;
641  bgLog("line\n");
642  break;
643  case FC_CUSTOM:
644  custx_ = lcustx_;
645  custy_ = lcusty_;
646  ncust_ = nlcust_;
647  fcomp = &BgEdgeDetect::CustomRegionComp;
648  feval = &BgEdgeDetect::CustomRegionEval;
649  bgLog("custom\n");
650  break;
651  }
652  confTr_ = (float) cl;
653  rankTr_ = (float) rl;
654  StrConfEstim(permNmxRank_, permNmxConf_, tdl, feval);
655 
656  //hysteresis thresholding
657 
658  grx_ = permGx_;
659  gry_ = permGy_;
660  int minpt = nMin;
661  NewHysteresisTr(tdh, tdl, cel, minpt, tr, tc);
662  bgLog("Done edge detection.\n");
663 
664  delete [] tdl;
665  delete [] tdh;
666  delete [] tr;
667  delete [] tc;
668 }
669 
670 void BgEdgeDetect::SaveNmxValues()
671 {
672  FILE* fd;
673  int i,j;
674 
675  fd = fopen("ranknmx.dat", "w");
676  for (j=0; j<y_; j++)
677  {
678  for (i=0; i<x_; i++)
679  {
680  fprintf(fd, "%f ", *(permNmxRank_+j*x_+i));
681  }
682  fprintf(fd, "\n");
683  }
684  fclose(fd);
685 
686  fd=fopen("confnmx.dat", "w");
687  for (j=0; j<y_; j++)
688  {
689  for (i=0; i<x_; i++)
690  {
691  fprintf(fd, "%f ", *(permNmxConf_+j*x_+i));
692  }
693  fprintf(fd, "\n");
694  }
695  fclose(fd);
696 }
697 
698 //Computes confedence map and rank
699 //Pre : cim is an image
700 //Post: confidence map and rank has been computed for cim
701 // and stored into confMap and rank respectively
702 void BgEdgeDetect::ComputeEdgeInfo(BgImage* cim, float* confMap, float *rank)
703 {
704  x_ = cim->x_;
705  y_ = cim->y_;
706  bgLog("Computing confidence map...\n");
707  float *pGx, *pGy, *pTemp;
708  BgImage tcim(x_, y_);
709  if (cim->colorIm_==1)
710  {
711  tcim.SetImageFromRGB(cim->im_, x_, y_, false);
712  } else
713  {
714  tcim = *cim;
715  }
716 
717  pGx = new float[x_*y_];
718  pGy = new float[x_*y_];
719  pTemp = new float[x_*y_];
720 
721  // compute gradient images
722  bgLog("...smooth-differentiation filtering\n");
723  GaussDiffFilter(&tcim, pGx, pGy, pTemp);
724 
725  // compute confidences (subspace estimate)
726  bgLog("...subspace estimate\n");
727  SubspaceEstim(pTemp, pGx, pGy, confMap);
728 
729  // compute edge strength from gradient image
730  bgLog("...edge strengths\n");
731  Strength(pGx, pGy, pTemp);
732 
733  // compute ranks of the strengths
734  bgLog("...computing ranks\n");
735  CompRanks(pTemp, rank);
736 
737  //de-allocate memory
738  delete [] pTemp;
739  delete [] pGy;
740  delete [] pGx;
741 }
742 /*
743 void BgEdgeDetect::ComputeConfidenceMap1(BgImage* cim, float* confMap)
744 {
745  ComputeConfidenceMap(cim, confMap);
746  BgEdgeList el;
747  DoEdgeDetect(cim, &el, RANK_NMX, CONF_NMX, RANK_H, CONF_H, RANK_L, CONF_L,
748  NMIN, FC_ELLIPSE, FC_SQUARE_BOX, FC_ELLIPSE);
749  BgImage tempImage(cim->x_, cim->y_);
750  el.SetBinImage(&tempImage);
751  int i;
752  for (i=0; i<(cim->x_*cim->y_); i++)
753  if (tempImage.im_[i] == 0)
754  confMap[i] = 0;
755 }
756 */
757 void BgEdgeDetect::DoEdgeDetect(BgImage* cim, BgEdgeList* cel, double nmxr, double nmxc,
758  double rh, double ch, double rl, double cl,
759  int nMin, int nmxType, int hystTypeHigh, int hystTypeLow)
760 {
761  x_ = cim->x_;
762  y_ = cim->y_;
763  bgLog("Start edge detection...\n");
764  permGx_ = new float[x_*y_];
765  permGy_ = new float[x_*y_];
766  permConf_ = new float[x_*y_];
767  permRank_ = new float[x_*y_];
768  permNmxRank_ = new float[x_*y_];
769  permNmxConf_ = new float[x_*y_];
770  havePerm_ = true;
771  float* tr;
772  float* tc;
773  float* tdh;
774  float* tdl;
775 
776  tr = new float[x_*y_];
777  tc = new float[x_*y_];
778  tdh = new float[x_*y_];
779  tdl = new float[x_*y_];
780 
781  // compute gradient images
782  bgLog("...smooth-differentiation filtering\n");
783  GaussDiffFilter(cim, permGx_, permGy_, tr);
784 
785  // compute confidences (subspace estimate)
786  bgLog("...subspace estimate\n");
787  SubspaceEstim(tr, permGx_, permGy_, permConf_);
788 
789  // compute edge strength from gradient image
790  bgLog("...edge strengths\n");
791  Strength(permGx_, permGy_, tr);
792 
793  // compute ranks of the strengths
794  bgLog("...computing ranks\n");
795  CompRanks(tr, permRank_);
796 
797  // new nonmaxima supression
798  bgLog("...nonmaxima supression: ");
799 
800  // select appropriate function
801  float (BgEdgeDetect::*fcomp)(float,float,float,float);
802  float (BgEdgeDetect::*feval)(float,float);
803  switch(nmxType)
804  {
805  case FC_ELLIPSE:
806  fcomp = &BgEdgeDetect::EllipseComp;
807  feval = &BgEdgeDetect::EllipseEval;
808  bgLog("arc\n");
809  break;
810  case FC_VERT_LINE:
811  fcomp = &BgEdgeDetect::VerticalLineComp;
812  feval = &BgEdgeDetect::VerticalLineEval;
813  bgLog("vertical line\n");
814  break;
815  case FC_HORIZ_LINE:
816  fcomp = &BgEdgeDetect::HorizontalLineComp;
817  feval = &BgEdgeDetect::HorizontalLineEval;
818  bgLog("horizontal line\n");
819  break;
820  case FC_SQUARE_BOX:
821  fcomp = &BgEdgeDetect::SquareComp;
822  feval = &BgEdgeDetect::SquareEval;
823  bgLog("box\n");
824  break;
825  case FC_LINE:
826  fcomp = &BgEdgeDetect::LineComp;
827  feval = &BgEdgeDetect::LineEval;
828  bgLog("line\n");
829  break;
830  default:
831  bgLog("Type not known\n");
832  return;
833  }
834 
835  confTr_ = (float) nmxc;
836  rankTr_ = (float) nmxr;
837  NewNonMaxSupress(permRank_, permConf_, permGx_, permGy_, permNmxRank_, permNmxConf_, fcomp);
838 
839  // new hysteresis thresholding
840  bgLog("...hysteresis thresholding, high: ");
841 
842  // select function, high curve
843  switch(hystTypeHigh)
844  {
845  case FC_ELLIPSE:
846  fcomp = &BgEdgeDetect::EllipseComp;
847  feval = &BgEdgeDetect::EllipseEval;
848  bgLog("arc");
849  break;
850  case FC_VERT_LINE:
851  fcomp = &BgEdgeDetect::VerticalLineComp;
852  feval = &BgEdgeDetect::VerticalLineEval;
853  bgLog("vertical line");
854  break;
855  case FC_HORIZ_LINE:
856  fcomp = &BgEdgeDetect::HorizontalLineComp;
857  feval = &BgEdgeDetect::HorizontalLineEval;
858  bgLog("horizontal line");
859  break;
860  case FC_SQUARE_BOX:
861  fcomp = &BgEdgeDetect::SquareComp;
862  feval = &BgEdgeDetect::SquareEval;
863  bgLog("box");
864  break;
865  case FC_LINE:
866  fcomp = &BgEdgeDetect::LineComp;
867  feval = &BgEdgeDetect::LineEval;
868  bgLog("line");
869  break;
870  case FC_CUSTOM:
871  custx_ = hcustx_;
872  custy_ = hcusty_;
873  ncust_ = nhcust_;
874  fcomp = &BgEdgeDetect::CustomRegionComp;
875  feval = &BgEdgeDetect::CustomRegionEval;
876  bgLog("custom");
877  break;
878  }
879 
880  confTr_ = (float) ch;
881  rankTr_ = (float) rh;
882  StrConfEstim(permNmxRank_, permNmxConf_, tdh, feval);
883 
884  bgLog(" low: ");
885 
886  // select function, low curve
887  switch(hystTypeLow)
888  {
889  case FC_ELLIPSE:
890  fcomp = &BgEdgeDetect::EllipseComp;
891  feval = &BgEdgeDetect::EllipseEval;
892  bgLog("arc\n");
893  break;
894  case FC_VERT_LINE:
895  fcomp = &BgEdgeDetect::VerticalLineComp;
896  feval = &BgEdgeDetect::VerticalLineEval;
897  bgLog("vertical line\n");
898  break;
899  case FC_HORIZ_LINE:
900  fcomp = &BgEdgeDetect::HorizontalLineComp;
901  feval = &BgEdgeDetect::HorizontalLineEval;
902  bgLog("horizontal line\n");
903  break;
904  case FC_SQUARE_BOX:
905  fcomp = &BgEdgeDetect::SquareComp;
906  feval = &BgEdgeDetect::SquareEval;
907  bgLog("box\n");
908  break;
909  case FC_LINE:
910  fcomp = &BgEdgeDetect::LineComp;
911  feval = &BgEdgeDetect::LineEval;
912  bgLog("line\n");
913  break;
914  case FC_CUSTOM:
915  custx_ = lcustx_;
916  custy_ = lcusty_;
917  ncust_ = nlcust_;
918  fcomp = &BgEdgeDetect::CustomRegionComp;
919  feval = &BgEdgeDetect::CustomRegionEval;
920  bgLog("custom\n");
921  break;
922  }
923  confTr_ = (float) cl;
924  rankTr_ = (float) rl;
925 
926  StrConfEstim(permNmxRank_, permNmxConf_, tdl, feval);
927 
928  grx_ = permGx_;
929  gry_ = permGy_;
930 
931  NewHysteresisTr(tdh, tdl, cel, nMin, tr, tc);
932 
933  bgLog("Done edge detection.\n");
934 
935  delete [] tdl;
936  delete [] tdh;
937  delete [] tr;
938  delete [] tc;
939 }
940 
941 void BgEdgeDetect::SubspaceEstim(float* im, float* grx, float* gry, float* cee)
942 {
943  // im original image
944  // grx, gry gradient of image
945  // cee confidence edge estimate
946 
947  float* itim;
948  float* itgx;
949  float* itgy;
950  float* itcee;
951  itim = im;
952  itgx = grx;
953  itgy = gry;
954  itcee = cee;
955 
956  double* tae;
957  double* ti;
958 
959  ti = new double[WW_*WW_];
960 
961  int i,j,l,c;
962  double v1;
963  double angleEdge;
964  int WW2 = WW_*WW_;
965 
966  itim += WL_*x_;
967  itgx += WL_*x_;
968  itgy += WL_*x_;
969  itcee += WL_*x_;
970 
971  for (j=WL_; j<(y_-WL_); j++)
972  {
973  for (i=0; i<WL_; i++)
974  itcee[i] = 0;
975  itim += WL_;
976  itgx += WL_;
977  itgy += WL_;
978  itcee += WL_;
979 
980 
981  for (i=WL_; i<(x_-WL_); i++, itim++, itgx++, itgy++, itcee++)
982  {
983  if ((fabs(*itgx)+fabs(*itgy))>TOL_E)
984  {
985  angleEdge = (-atan2(*itgx, *itgy))*180.0/PI;
986  tae = lookTable_[(int) (angleEdge+180.49)];
987 
988  //A=A-mean(A)
989  v1=0;
990  for (l=0; l<WW_; l++)
991  {
992  for (c=0; c<WW_; c++)
993  {
994  v1 += ti[l*WW_+c] = *(itim+(l-WL_)*x_+c-WL_);
995  }
996  }
997  v1 /= WW2;
998  for (l=0; l<WW2; l++)
999  ti[l] -= v1;
1000 
1001  //A/norm(A,'fro')
1002  v1 = 0;
1003  for (l=0; l<WW2; l++)
1004  v1 += ti[l]*ti[l];
1005  v1 = sqrt(v1);
1006  for (l=0; l<WW2; l++)
1007  ti[l] /= v1;
1008 
1009  //global
1010  v1 = 0;
1011  for (l=0; l<WW2; l++)
1012  v1 += tae[l]*ti[l];
1013  v1 = fabs(v1);
1014  *itcee = (float) v1;
1015  }
1016  else
1017  {
1018  *itcee = 0;
1019  }
1020  }
1021  for (i=0; i<WL_; i++)
1022  itcee[i] = 0;
1023  itim += WL_;
1024  itgx += WL_;
1025  itgy += WL_;
1026  itcee += WL_;
1027  }
1028  WW2 = x_*y_;
1029  for (j=0; j<(WL_*x_); j++)
1030  {
1031  cee[j] = 0;
1032  cee[WW2-j-1] = 0;
1033  }
1034 
1035 
1036  delete [] ti;
1037 }
1038 
1039 void BgEdgeDetect::GaussFilter(BgImage* cim, float* fim, double sigma, int width)
1040 {
1041  double* filter;
1042  unsigned char* im;
1043  double* tim;
1044  double sum = 0;
1045  double sum1 = 0;
1046  int i,ii,jj,j,k;
1047 
1048  im = cim->im_;
1049  if (width==-2)
1050  {
1051  for (i=0; i<x_*y_;i++)
1052  fim[i] = im[i];
1053  return;
1054  }
1055 
1056  if (width<3)
1057  width = (int) (1+2*ceil(2.5*sigma));
1058  int tail = width/2;
1059  width = 2*tail+1;
1060 
1061  //create kernel
1062  filter = new double[width];
1063  for (i=-tail; i<=tail; i++)
1064  sum += filter[i+tail] = exp(-i*i/(2*sigma*sigma));
1065  for (i=0; i<width; i++)
1066  filter[i] /= sum;
1067 
1068  //filter image
1069  im = cim->im_;
1070  tim = new double[x_*y_];
1071  for (j=0; j<y_; j++)
1072  {
1073  for (i=tail; i<(x_-tail); i++)
1074  {
1075  sum=0;
1076  for (k=-tail; k<=tail; k++)
1077  sum += im[j*x_+i+k]*filter[k+tail];
1078  tim[j*x_+i] = sum;
1079  }
1080 
1081  for (i=0; i<tail; i++)
1082  {
1083  tim[j*x_+i] = 0;
1084  tim[j*x_+x_-i-1] = 0;
1085  for (k=-tail; k<=tail; k++)
1086  {
1087  ii = (k+i)>=0 ? k+i : 0;
1088  tim[j*x_+i] += im[j*x_+ii]*filter[k+tail];
1089  ii = (x_-i-1+k)<x_ ? x_-i-1+k : x_-1;
1090  tim[j*x_+x_-i-1] += im[j*x_+ii]*filter[k+tail];
1091  }
1092  }
1093  }
1094 
1095  for (i=0; i<x_; i++)
1096  {
1097  for (j=tail; j<(y_-tail); j++)
1098  {
1099  sum=0;
1100  for (k=-tail; k<=tail; k++)
1101  sum += tim[(j+k)*x_+i]*filter[k+tail];
1102  fim[j*x_+i] = (float) (sum);
1103  }
1104  for (j=0; j<tail; j++)
1105  {
1106  sum = 0;
1107  sum1 = 0;
1108  for (k=-tail; k<=tail; k++)
1109  {
1110  jj = (k+j)>=0 ? k+j : 0;
1111  sum += tim[jj*x_+i]*filter[k+tail];
1112  jj = (y_-j-1+k)<y_ ? y_-j-1+k : y_-1;
1113  sum1 += tim[jj*x_+i]*filter[k+tail];
1114  }
1115  fim[j*x_+i] = (float) (sum);
1116  fim[(y_-j-1)*x_+i] = (float) (sum1);
1117  }
1118  }
1119  delete [] filter;
1120  delete [] tim;
1121 }
1122 
1123 void BgEdgeDetect::GaussDiffFilter(BgImage* cim, float* grx, float* gry, float* rezIm)
1124 {
1125 
1126  double* sf; //smooth filter
1127  double* df; //diff filter
1128  unsigned char* im;
1129 
1130  double* tim;
1131  double sum = 0;
1132  double sum1 = 0;
1133  int i, j, k;
1134 
1135  //create kernels
1136  sf = smofil_;
1137  df = diffil_;
1138 
1139  im = cim->im_;
1140  tim = new double[x_*y_];
1141  for (i=0; i<x_*y_; i++)
1142  {
1143  grx[i] = gry[i] = 0;
1144  tim[i] = 0;
1145  rezIm[i] = im[i];
1146  }
1147 
1148  //filter image x
1149  //smooth on y
1150  for (i=0; i<x_; i++)
1151  {
1152  for (j=WL_; j<(y_-WL_); j++)
1153  {
1154  sum = 0;
1155  for (k=-WL_; k<=WL_; k++)
1156  sum += im[(j+k)*x_+i]*sf[k+WL_];
1157  tim[j*x_+i] = sum;
1158  }
1159  }
1160  //diff on x
1161  for (j=0; j<y_; j++)
1162  {
1163  for (i=WL_; i<(x_-WL_); i++)
1164  {
1165  sum = 0;
1166  for (k=-WL_; k<=WL_; k++)
1167  sum += tim[j*x_+i+k]*df[k+WL_];
1168  grx[j*x_+i] = (float) (sum);
1169  }
1170  }
1171 
1172  //filter image y
1173  for (i=0; i<x_*y_;i++)
1174  tim[i] = 0;
1175  im = cim->im_;
1176  //smooth on x
1177  for (j=0; j<y_; j++)
1178  {
1179  for (i=WL_; i<(x_-WL_); i++)
1180  {
1181  sum = 0;
1182  for (k=-WL_; k<=WL_; k++)
1183  sum += im[j*x_+i+k]*sf[k+WL_];
1184  tim[j*x_+i] = sum;
1185  }
1186  }
1187  //diff on y
1188  for (i=0; i<x_; i++)
1189  {
1190  for (j=WL_; j<(y_-WL_); j++)
1191  {
1192  sum = 0;
1193  for (k=-WL_; k<=WL_; k++)
1194  sum += tim[(j+k)*x_+i]*df[k+WL_];
1195  gry[j*x_+i] = (float) (sum);
1196  }
1197  }
1198  delete [] tim;
1199 }
1200 
1201 void BgEdgeDetect::CompRanks(float* strength, float* ranks)
1202 {
1203  int* index;
1204  float* ra;
1205  ra = new float[x_*y_];
1206  index = new int[x_*y_];
1207  int ii;
1208 
1209  for (ii=0; ii<x_*y_; ii++)
1210  {
1211  index[ii] = ii;
1212  ranks[ii] = 0;
1213  ra[ii] = strength[ii];
1214  }
1215 
1216  //heap sort with ranks (from numerical recipies)
1217  unsigned long i, ir, j, l;
1218  unsigned long n;
1219  n = x_*y_;
1220  float rra;
1221  int irra;
1222 
1223  if (n<2)
1224  return;
1225  l = (n>>1)+1;
1226  ir = n;
1227  for (;;)
1228  {
1229  if (l>1)
1230  {
1231  rra = ra[--l-1];
1232  irra = index[l-1];
1233  }
1234  else
1235  {
1236  rra = ra[ir-1];
1237  irra = index[ir-1];
1238  ra[ir-1] = ra[1-1];
1239  index[ir-1] = index[1-1];
1240  if (--ir==1)
1241  {
1242  ra[1-1] = rra;
1243  index[1-1] = irra;
1244  break;
1245  }
1246  }
1247  i = l;
1248  j = l+l;
1249  while (j<=ir)
1250  {
1251  if (j<ir && ra[j-1]<ra[j+1-1])
1252  j++;
1253  if (rra<ra[j-1])
1254  {
1255  ra[i-1] = ra[j-1];
1256  index[i-1] = index[j-1];
1257  i = j;
1258  j <<= 1;
1259  }
1260  else
1261  j = ir+1;
1262  }
1263  ra[i-1] = rra;
1264  index[i-1] = irra;
1265  }
1266 
1267  //setranks
1268  irra = 1;
1269  for (ii=1; ii<x_*y_; ii++)
1270  {
1271  if (ra[ii]>ZERO_TRESH)
1272  {
1273  ranks[index[ii]] = (float) irra;
1274  if (ra[ii]>ra[ii-1])
1275  irra++;
1276  }
1277  }
1278  irra--;
1279  for (ii=0; ii<x_*y_; ii++)
1280  ranks[ii] /= irra;
1281 
1282  delete [] index;
1283  delete [] ra;
1284 }
1285 
1286 void BgEdgeDetect::StrConfEstim(float* ranks, float* confidence, float* rezult,
1287  float (BgEdgeDetect::*feval)(float,float))
1288 {
1289  int i;
1290  for (i=0; i<x_*y_; i++)
1291  {
1292  rezult[i] = (this->*feval)(ranks[i], confidence[i]);
1293  }
1294 }
1295 
1296 void BgEdgeDetect::Strength(float* grx, float* gry, float* strength)
1297 {
1298  int i,j;
1299  float* itgx;
1300  float* itgy;
1301  float* its;
1302  double val;
1303  itgx = grx;
1304  itgy = gry;
1305  its = strength;
1306  for (j=0; j<y_; j++)
1307  for(i=0; i<x_; i++)
1308  {
1309  val = sqrt(((double) (*itgx*(*(itgx++))))+((double) (*itgy*(*(itgy++)))));
1310  *(its++)=(float) (val);
1311  }
1312 }
1313 
1314 void BgEdgeDetect::NewNonMaxSupress(float* rank, float* conf, float* grx, float* gry, float* nmxRank, float* nmxConf,
1315  float (BgEdgeDetect::*feval)(float, float, float, float))
1316 {
1317  int i,j;
1318  float* itr;
1319  float* itc;
1320  float* itgx;
1321  float* itgy;
1322  float* itnmxr;
1323  float* itnmxc;
1324  float alpha,r1,c1,r2,c2,lambda;
1325  itr = rank;
1326  itc = conf;
1327  itgx = grx;
1328  itgy = gry;
1329  itnmxr = nmxRank;
1330  itnmxc = nmxConf;
1331 
1332  for (i=0; i<x_*y_; i++)
1333  {
1334  itnmxr[i] = itnmxc[i] = 0;
1335  }
1336  for(i=0; i<x_; i++)
1337  {
1338  itr[i] = itc[i] = 0;
1339  itr[(y_-1)*x_+i] = itc[(y_-1)*x_+i] = 0;
1340  }
1341  for(j=0; j<y_; j++)
1342  {
1343  itr[j*x_] = itc[j*x_] = 0;
1344  itr[j*x_+x_-1] = itc[j*x_+x_-1] = 0;
1345  }
1346 
1347  for (j=0; j<y_; j++)
1348  {
1349  for (i=0; i<x_; i++, itr++, itc++, itgx++, itgy++, itnmxr++, itnmxc++)
1350  {
1351  if (*itr>0 && *itc>0)
1352  {
1353  alpha = (float) atan2(*itgy, *itgx);
1354  alpha = (alpha<0) ? alpha+(float)PI : alpha;
1355  if (alpha<=PI/4)
1356  {
1357  lambda = (float) tan(alpha);
1358  r1 = (1-lambda)*(*(itr+1))+lambda*(*(itr+x_+1));
1359  c1 = (1-lambda)*(*(itc+1))+lambda*(*(itc+x_+1));
1360  r2 = (1-lambda)*(*(itr-1))+lambda*(*(itr-x_-1));
1361  c2 = (1-lambda)*(*(itc-1))+lambda*(*(itc-x_-1));
1362  if ((this->*feval)(*itr, *itc, r1, c1)<0 && (this->*feval)(*itr, *itc, r2, c2)<=0)
1363  {
1364  *itnmxr = *itr;
1365  *itnmxc = *itc;
1366  }
1367  }
1368  else if (alpha<=PI/2)
1369  {
1370  lambda = (float) tan(PI/2-alpha);
1371  r1 = (1-lambda)*(*(itr+x_))+lambda*(*(itr+x_+1));
1372  c1 = (1-lambda)*(*(itc+x_))+lambda*(*(itc+x_+1));
1373  r2 = (1-lambda)*(*(itr-x_))+lambda*(*(itr-x_-1));
1374  c2 = (1-lambda)*(*(itc-x_))+lambda*(*(itc-x_-1));
1375  if ((this->*feval)(*itr, *itc, r1, c1)<0 && (this->*feval)(*itr, *itc, r2, c2)<=0)
1376  {
1377  *itnmxr = *itr;
1378  *itnmxc = *itc;
1379  }
1380 
1381  }
1382  else if(alpha<=3*PI/4)
1383  {
1384  lambda = (float) tan(alpha-PI/2);
1385  r1 = (1-lambda)*(*(itr+x_))+lambda*(*(itr+x_-1));
1386  c1 = (1-lambda)*(*(itc+x_))+lambda*(*(itc+x_-1));
1387  r2 = (1-lambda)*(*(itr-x_))+lambda*(*(itr-x_+1));
1388  c2 = (1-lambda)*(*(itc-x_))+lambda*(*(itc-x_+1));
1389  if ((this->*feval)(*itr, *itc, r1, c1)<0 && (this->*feval)(*itr, *itc, r2, c2)<=0)
1390  {
1391  *itnmxr = *itr;
1392  *itnmxc = *itc;
1393  }
1394  }
1395  else
1396  {
1397  lambda = (float) tan(PI-alpha);
1398  r1 = (1-lambda)*(*(itr-1))+lambda*(*(itr+x_-1));
1399  c1 = (1-lambda)*(*(itc-1))+lambda*(*(itc+x_-1));
1400  r2 = (1-lambda)*(*(itr+1))+lambda*(*(itr-x_+1));
1401  c2 = (1-lambda)*(*(itc+1))+lambda*(*(itc-x_+1));
1402  if ((this->*feval)(*itr, *itc, r1, c1)<0 && (this->*feval)(*itr, *itc, r2, c2)<=0)
1403  {
1404  *itnmxr = *itr;
1405  *itnmxc = *itc;
1406  }
1407  }
1408  }
1409  }
1410  }
1411 }
1412 
1413 void BgEdgeDetect::NewHysteresisTr(float* edge, float* low, BgEdgeList* cel, int nMin, float* mark, float* coord)
1414 {
1415  float* tm;
1416  float* te;
1417  float* tl;
1418  int i,j,n;
1419 
1420  n=0;
1421  for (i=0, tm=mark; i<x_*y_; i++,tm++)
1422  *tm=0;
1423 
1424  te_ = te = edge;
1425  tm_ = tm = mark;
1426  tl_ = tl = low;
1427 
1428  for (j=0; j<y_; j++)
1429  {
1430  for (i=0; i<x_; i++, tm++, te++)
1431  {
1432  if ((*tm==0) && ((*te)>HYST_LOW_CUT))
1433  {
1434  //found an edge start
1435  npt_ = 0;
1436  *tm = 1;
1437  tc_ = coord;
1438  NewEdgeFollow(i, j);
1439  //store the edge
1440  if (npt_>=nMin) cel->AddEdge(coord, npt_);
1441  }
1442  }
1443  }
1444 }
1445 
1446 void BgEdgeDetect::NewEdgeFollow(int ii,int jj)
1447 {
1448  int i;
1449  int iin, jjn;
1450  for (i=0; i<8; i++)
1451  {
1452  iin = ii+gNb[i][0];
1453  jjn = jj+gNb[i][1];
1454  if ((tm_[jjn*x_+iin]==0) && ((tl_[jjn*x_+iin])>0))
1455  {
1456  tm_[jjn*x_+iin] = 1;
1457  NewEdgeFollow(iin, jjn);
1458  }
1459  }
1460  *(tc_++) = (float) ii;
1461  *(tc_++) = (float) jj;
1462  npt_++;
1463 }
1464 
1465 void BgEdgeDetect::SetCustomHigh(int* x, int* y, int n, int sx, int sy)
1466 {
1467  if (nhcust_>0)
1468  {
1469  delete [] hcustx_;
1470  delete [] hcusty_;
1471  }
1472  nhcust_ = 0;
1473  hcustx_ = hcusty_ = 0;
1474  nhcust_ = n+2;
1475  hcustx_ = new float[nhcust_];
1476  hcusty_ = new float[nhcust_];
1477 
1478  int idx,i;
1479  idx = 0;
1480  hcustx_[idx] = 0;
1481  hcusty_[idx++] = 0;
1482  for (i=0; i<n; i++)
1483  {
1484  hcustx_[idx] = ((float) x[i])/sx;
1485  hcusty_[idx++] = (float)(1.0-((float) y[i])/sy);
1486  }
1487  hcustx_[idx] = 0;
1488  hcusty_[idx++] = 0;
1489 
1490  bgLog(" hyst high custom x: ");
1491  for (i=0; i<=n; i++)
1492  bgLog(" %f", hcustx_[i]);
1493  bgLog("\n");
1494  bgLog(" hist high custom y: ");
1495  for (i=0; i<=n; i++)
1496  bgLog(" %f", hcusty_[i]);
1497  bgLog("\n");
1498 }
1499 
1500 void BgEdgeDetect::SetCustomHigh(double* x, double* y, int n)
1501 {
1502  if (nhcust_>0)
1503  {
1504  delete [] hcustx_;
1505  delete [] hcusty_;
1506  }
1507  nhcust_ = 0;
1508  hcustx_ = hcusty_ = 0;
1509  nhcust_ = n+2;
1510  hcustx_ = new float[nhcust_];
1511  hcusty_ = new float[nhcust_];
1512 
1513  int idx,i;
1514  idx = 0;
1515  hcustx_[idx] = 0;
1516  hcusty_[idx++] = 0;
1517  for (i=0; i<n; i++)
1518  {
1519  hcustx_[idx] = (float) x[i];
1520  hcusty_[idx++] = (float) y[i];
1521  }
1522  hcustx_[idx] = 0;
1523  hcusty_[idx++] = 0;
1524 
1525  bgLog(" hyst high custom x: ");
1526  for (i=0; i<=n; i++)
1527  bgLog(" %f", hcustx_[i]);
1528  bgLog("\n");
1529  bgLog(" hist high custom y: ");
1530  for (i=0; i<=n; i++)
1531  bgLog(" %f", hcusty_[i]);
1532  bgLog("\n");
1533 }
1534 
1535 void BgEdgeDetect::SetCustomLow(int* x, int* y, int n, int sx, int sy)
1536 {
1537  if(nlcust_>0)
1538  {
1539  delete [] lcustx_;
1540  delete [] lcusty_;
1541  }
1542  nlcust_ = 0;
1543  lcustx_ = lcusty_ = 0;
1544  nlcust_ = n+2;
1545  lcustx_ = new float[nlcust_];
1546  lcusty_ = new float[nlcust_];
1547 
1548  int idx,i;
1549  idx = 0;
1550  lcustx_[idx] = 0;
1551  lcusty_[idx++] = 0;
1552  for (i=0; i<n; i++)
1553  {
1554  lcustx_[idx] = ((float) x[i])/sx;
1555  lcusty_[idx++] = (float)(1.0-((float) y[i])/sy);
1556  }
1557  lcustx_[idx] = 0;
1558  lcusty_[idx++] = 0;
1559  bgLog(" hyst low custom x: ");
1560  for (i=0; i<=n; i++)
1561  bgLog(" %f", lcustx_[i]);
1562  bgLog("\n");
1563  bgLog(" low custom y: ");
1564  for (i=0; i<=n; i++)
1565  bgLog(" %f", lcusty_[i]);
1566  bgLog("\n");
1567 }
1568 
1569 void BgEdgeDetect::SetCustomLow(double* x, double* y, int n)
1570 {
1571  if(nlcust_>0)
1572  {
1573  delete [] lcustx_;
1574  delete [] lcusty_;
1575  }
1576  nlcust_ = 0;
1577  lcustx_ = lcusty_ = 0;
1578  nlcust_ = n+2;
1579  lcustx_ = new float[nlcust_];
1580  lcusty_ = new float[nlcust_];
1581 
1582  int idx,i;
1583  idx = 0;
1584  lcustx_[idx] = 0;
1585  lcusty_[idx++] = 0;
1586  for (i=0; i<n; i++)
1587  {
1588  lcustx_[idx] = (float) x[i];
1589  lcusty_[idx++] = (float) y[i];
1590  }
1591  lcustx_[idx] = 0;
1592  lcusty_[idx++] = 0;
1593  bgLog(" hyst low custom x: ");
1594  for (i=0; i<=n; i++)
1595  bgLog(" %f", lcustx_[i]);
1596  bgLog("\n");
1597  bgLog(" low custom y: ");
1598  for (i=0; i<=n; i++)
1599  bgLog(" %f", lcusty_[i]);
1600  bgLog("\n");
1601 }