24 #include "../ami_image/image.h"
25 #include "../ami_utilities/utilities.h"
26 #include "../ami_primitives/subpixel_image_contours.h"
36 template <
class T,
class U>
46 if(img.get_roi().size()==6){
47 img.get_roi_image(img_conv);
50 if(img.width()!=img_conv.width() || img.height()!=img_conv.height() ||
51 img.nChannels()!=img_conv.nChannels()){
52 img_conv=
ami::image<U>(img.width(),img.height(),img.nChannels());
54 int k,k_end=img.width()*img.height()*img.nChannels();
56 #pragma omp parallel for shared(k_end) private(k)
58 for(k=0;k<k_end;k++) img_conv[k]= (U) img[k];
61 unsigned int width=img_conv.width();
62 unsigned int height=img_conv.height();
63 unsigned int nChannels=img_conv.nChannels();
64 unsigned int image_size=width*height;
65 unsigned int image_size_total=image_size*nChannels;
67 if(precision<0) precision=0;
69 unsigned int Nc_x=(
unsigned int) (precision*sigma_x>1?precision*sigma_x:1);
70 unsigned int Nc_y=(
unsigned int) (precision*sigma_y>1?precision*sigma_y:1);
71 unsigned int Nc_a=(Nc_x>Nc_y)?Nc_x:Nc_y;
73 double t_y=sigma_x*sigma_x/(2*Nc_x);
74 double t_x=sigma_y*sigma_y/(2*Nc_y);
75 double l_x,v_x,l_y,v_y,l_x_1,l_y_1;
78 if(t_x==0 && t_y==0)
return;
82 l_x=(float)(1.+2.*t_x-sqrt((
double) 4*t_x+1))/(2*t_x);
87 l_y=(float)(1.+2.*t_y-sqrt((
double) 4*t_y+1))/(2*t_y);
94 int c,y,x,x_end,m,y_end;
95 for(cont=0;cont<(int)Nc_a;cont++){
98 #pragma omp parallel \
99 shared(width,image_size_total,image_size,l_y,v_y,l_y_1) \
100 private(c,x,y,m,x_end)
102 for(c=0;c<(int)image_size_total;c+=image_size){
104 #pragma omp for nowait
106 for(y=0;y<(int)image_size;y+=width){
108 vector <U> paso(width);
109 paso[0]= img_conv[m]/l_y_1;
112 for(x=m+1;x<x_end;x++,n++){
113 paso[n]= img_conv[x]+l_y*paso[n-1];
115 img_conv[x_end-1]=paso[width-1]/l_y_1;
117 for(x=x_end-2;x>m;x--,n--){
118 img_conv[x]= paso[n]+l_y*img_conv[x+1];
120 img_conv[x]= paso[n]+l_y*img_conv[x+1];
121 for(x=m;x<x_end;x++){
129 #pragma omp parallel \
130 shared(width,height,image_size_total,image_size,l_x,v_x,l_x_1) \
131 private(c,x,y,m,y_end)
133 for(c=0;c<(int)image_size_total;c+=image_size){
135 #pragma omp for nowait
137 for(x=0;x<(int)width;x++){
139 vector <U> paso(height);
140 paso[0]= img_conv[m]/l_x_1;
143 for(y=m+width;y<y_end;y+=width,n++){
144 paso[n]= img_conv[y]+l_x*paso[n-1];
146 img_conv[y_end-width]=paso[height-1]/l_x_1;
148 for(y=y_end-2*width;y>m;y-=width,n--){
149 img_conv[y]= paso[n]+l_x*img_conv[y+width];
151 img_conv[y]= paso[n]+l_x*img_conv[y+width];
152 for(y=m;y<y_end;y+=width){
169 template <
class T,
class U>
174 const bool NeigborhoodType )
177 if(img.get_roi().size()==6){
179 img.get_roi_image(img2);
180 grad(img2,grad_x,grad_y,NeigborhoodType);
185 unsigned int width=img.width();
186 unsigned int height=img.height();
187 unsigned int nChannels=img.nChannels();
188 unsigned int size_image=width*height;
189 unsigned int size_image_width=size_image-width;
190 unsigned int total_image_size=size_image*nChannels;
191 unsigned int width_1=width-1;
193 if((
int)width!=grad_x.width() || (int)height!=grad_x.height() ||
194 (int)nChannels!=grad_x.nChannels()){
198 if((
int)width!=grad_y.width() || (int)height!=grad_y.height() ||
199 (int)nChannels!=grad_y.nChannels()){
203 if(NeigborhoodType==0){
206 vector < vector <unsigned int> > b=boundary_neighborhood_5n(width,height);
207 int m,c,y,k,k_end,b_size=b.size();
210 #pragma omp parallel \
211 shared(width,width_1,total_image_size,size_image,size_image_width,b,b_size)\
212 private(c,y,m,k,k_end)
214 for(c=0;c<(int)total_image_size;c+=size_image){
216 #pragma omp for nowait
218 for(y=width;y<(int)size_image_width;y+=width){
221 for(k=m+1; k<k_end; k++){
222 grad_y[k]=(U) img[k+width]- (U) img[k-width];
223 grad_x[k]=(U) img[k-1]- (U) img[k+1];
227 #pragma omp for nowait
229 for(
int k=0;k<b_size;k++){
230 grad_y[b[k][0]+c]=(U) img[b[k][2]+c]- (U) img[b[k][1]+c];
231 grad_x[b[k][0]+c]=(U) img[b[k][4]+c]- (U) img[b[k][3]+c];
236 vector < vector <unsigned int> > b=boundary_neighborhood_9n(width,height);
237 double coef1,coef2,c1,d1;
238 coef1=sqrt((
double) 2.);
239 coef2=0.25*(2.-coef1);
241 int m,c,y,k,l,k_end,b_size=b.size();
243 #pragma omp parallel \
244 shared(width,width_1,total_image_size,size_image,size_image_width,b,b_size,\
245 coef1,coef2) private(c,y,m,k,k_end,c1,d1)
247 for(c=0;c<(int)total_image_size;c+=size_image){
249 #pragma omp for nowait
251 for(y=width;y<(int)size_image_width;y+=width){
254 for(k=m+1; k<k_end; k++){
255 c1=img[k+width+1]-img[k-width-1];
256 d1=img[k-width+1]-img[k+width-1];
257 grad_y[k]=(U)(coef1*((U) img[k+width]-(U)img[k-width])+coef2*(c1-d1));
258 grad_x[k]=(U)(-(coef1*((U) img[k+1]- (U) img[k-1])+coef2*(c1+d1)));
262 #pragma omp for nowait
264 for(l=0;l<b_size;l++){
265 c1=img[c+b[l][6]]-img[c+b[l][7]];
266 d1=img[c+b[l][5]]-img[c+b[l][8]];
267 grad_y[c+b[l][0]]=(U)(coef1*((U) img[c+b[l][2]]- (U) img[c+b[l][1]])+
269 grad_x[c+b[l][0]]=(U)(-(coef1*((U) img[c+b[l][3]]- (U) img[c+b[l][4]])+
289 ami_malloc1d(y,
float,n);
290 for(
int mm=0;mm<n;mm++) y[mm]=x[mm];
296 if (ir == l+1 && y[ir] < y[l]) {
297 paso=y[l]; y[l]=y[ir]; y[ir]=paso;
304 paso=y[mid]; y[mid]=y[l+1]; y[l+1]=paso;
306 paso=y[l]; y[l]=y[ir]; y[ir]=paso;
308 if (y[l+1] > y[ir]) {
309 paso=y[l+1]; y[l+1]=y[ir]; y[ir]=paso;
312 paso=y[l]; y[l]=y[l+1]; y[l+1]=paso;
319 do i++;
while (y[i] < a);
320 do j--;
while (y[j] > a);
322 paso=y[i]; y[i]=y[j]; y[j]=paso;
336 double low_threshold,
339 if(imap[k]==1)
return;
340 int width = imap.width();
343 int nv[8] = {k-1, k+1, k-width, k+width,
344 k-width-1, k-width+1, k+width-1, k+width+1};
346 for(
int iv=0; iv<8; iv++)
349 if(imap[nk]==3)
continue;
350 if((imap[nk]==0) && (gradien_norm[nk]>low_threshold))
351 fill_imap(imap, gradien_norm, low_threshold, nk);
375 bool neighborhoodtype = 9;
376 int width = input.width();
377 int height = input.height();
378 int size = width*height;
380 float low_threshold,high_threshold;
384 float *gradient_norm =
new float[size];
395 grad(blurred,x_grad,y_grad,neighborhoodtype);
396 for(
int k=0; k<size; k++)
397 gradient_norm[k] = sqrt((x_grad[k]*x_grad[k]) + (y_grad[k]*y_grad[k]));
403 low_threshold=(low_threshold<2)?2:low_threshold;
404 high_threshold=(high_threshold<4)?4:high_threshold;
407 for(
int i=0; i<height; i++)
409 for(
int j=0; j<width; j++)
413 if(gradient_norm[pos] > high_threshold)
418 if((i==0) || (i==height-1) || (j==0) || (j==width-1))
426 float c = 1/(sqrt(2.));
427 for(
int i=1; i<height-1; i++)
429 for(
int j=1; j<width-1; j++)
433 float o1 = fabs(x_grad[k]);
434 float o2 = fabs(y_grad[k]);
435 float o3 = fabs(c*(x_grad[k]+y_grad[k]));
436 float o4 = fabs(c*(x_grad[k]-y_grad[k]));
441 float candidate = max(o1,max(o2,max(o3,o4)));
447 neigh1 = gradient_norm[k-1];
448 neigh2 = gradient_norm[k+1];
454 neigh1 = gradient_norm[k-width];
455 neigh2 = gradient_norm[k+width];
461 neigh1 = gradient_norm[k-width+1];
462 neigh2 = gradient_norm[k+width-1];
468 neigh1 = gradient_norm[k-width-1];
469 neigh2 = gradient_norm[k+width+1];
473 float p_grad = gradient_norm[k];
474 if((p_grad>neigh1) && (p_grad>neigh2))
476 if(gradient_norm[k] >= high_threshold)
489 for(
int k=0; k<size; k++)
494 int nv[8] = {k-1, k+1, k-width, k+width,
495 k-width-1, k-width+1, k+width-1, k+width+1};
496 for(
int iv=0; iv<8; iv++)
499 if(imap[nk]>1)
continue;
500 if((imap[nk]==0) && (gradient_norm[nk]>low_threshold))
502 fill_imap(imap, gradient_norm, low_threshold, nk);
518 for(
int i=1; i<height-1; i++)
520 for(
int j=1; j<width-1; j++)
535 coseno[pos] = x_grad[pos]/gradient_norm[pos];
536 seno[pos] = y_grad[pos]/gradient_norm[pos];
551 delete []gradient_norm;
565 float canny_low_threshold ,
566 float canny_high_threshold )
568 int size_=input.width()*input.height();
569 float *seno =
new float[size_];
570 float *coseno =
new float[size_];
571 int *x_pos =
new int[size_];
572 int *y_pos =
new int[size_];
574 canny(input,edges,seno,coseno,x_pos,y_pos,canny_low_threshold,canny_high_threshold);
578 for(
int i=0; i<size_; i++){
580 contours.
get_c()[i] =
true;
581 contours.get_x()[i] = (float)x_pos[i];
582 contours.get_y()[i] = (float)y_pos[i];
583 contours.get_coseno()[i] = seno[i];
584 contours.get_seno()[i] = coseno[i];
587 contours.get_c()[i] =
false;
Class to store subpixel contours.
Definition: subpixel_image_contours.h:39
void grad(const ami::image< T > &img, ami::image< U > &grad_x, ami::image< U > &grad_y, const bool NeigborhoodType)
Definition: filters.h:170
void canny(ami::image< T > input, ami::image< T > &output, float *seno, float *coseno, int *x, int *y, float per_low, float per_high)
Definition: filters.h:365
bool * get_c()
Return array c to identity contour points.
Definition: subpixel_image_contours.h:106
void gauss_conv(ami::image< T > &img, ami::image< U > &img_conv, double sigma_x, double sigma_y, double precision)
Definition: filters.h:37
float ami_median_float(int k, int n, float *x)
FUNCTION TO COMPUTE THE MEDIAN OF A VECTOR IN FLOAT PRECISION.
Definition: filters.h:285
Class to store multiChannel images and basic methods.
Definition: image.h:65