ElSe.cpp 38.4 KB
Newer Older
Thiago Santini's avatar
Thiago Santini committed
1 2 3 4
#include "ElSe.h"

#include <QDebug>
#include <QThread>
5
#include <opencv2/highgui.hpp>
Thiago Santini's avatar
Thiago Santini committed
6 7 8 9

using namespace cv;

std::string ElSe::desc = "ElSe (Fuhl et al. 2016)";
10 11
float ElSe::minArea = 0;
float ElSe::maxArea = 0;
Thiago Santini's avatar
Thiago Santini committed
12 13 14

#define IMG_SIZE 680 //400

15 16
static bool is_good_ellipse_eval(RotatedRect* ellipse, Mat* pic, int* erg)
{
Thiago Santini's avatar
Thiago Santini committed
17

18
    if (ellipse->center.x == 0 && ellipse->center.y == 0)
Thiago Santini's avatar
Thiago Santini committed
19 20
        return false;

21 22
    float x0 = ellipse->center.x;
    float y0 = ellipse->center.y;
Thiago Santini's avatar
Thiago Santini committed
23

24 25 26 27
    int st_x = (int)ceil(x0 - (ellipse->size.width / 4.0));
    int st_y = (int)ceil(y0 - (ellipse->size.height / 4.0));
    int en_x = (int)floor(x0 + (ellipse->size.width / 4.0));
    int en_y = (int)floor(y0 + (ellipse->size.height / 4.0));
Thiago Santini's avatar
Thiago Santini committed
28

29 30 31
    float val = 0.0;
    float val_cnt = 0;
    float ext_val = 0.0;
Thiago Santini's avatar
Thiago Santini committed
32

33 34
    for (int i = st_x; i < en_x; i++)
        for (int j = st_y; j < en_y; j++) {
Thiago Santini's avatar
Thiago Santini committed
35

36 37
            if (i > 0 && i < pic->cols && j > 0 && j < pic->rows) {
                val += pic->data[(pic->cols * j) + i];
Thiago Santini's avatar
Thiago Santini committed
38 39 40 41
                val_cnt++;
            }
        }

42 43
    if (val_cnt > 0)
        val = val / val_cnt;
Thiago Santini's avatar
Thiago Santini committed
44 45 46
    else
        return false;

47
    val_cnt = 0;
Thiago Santini's avatar
Thiago Santini committed
48

49 50 51 52
    st_x = (int)(x0 - (ellipse->size.width * 0.75));
    st_y = (int)(y0 - (ellipse->size.height * 0.75));
    en_x = (int)(x0 + (ellipse->size.width * 0.75));
    en_y = (int)(y0 + (ellipse->size.height * 0.75));
Thiago Santini's avatar
Thiago Santini committed
53

54 55 56 57
    int in_st_x = (int)ceil(x0 - (ellipse->size.width / 2));
    int in_st_y = (int)ceil(y0 - (ellipse->size.height / 2));
    int in_en_x = (int)floor(x0 + (ellipse->size.width / 2));
    int in_en_y = (int)floor(y0 + (ellipse->size.height / 2));
Thiago Santini's avatar
Thiago Santini committed
58

59 60 61 62 63 64 65
    for (int i = st_x; i < en_x; i++)
        for (int j = st_y; j < en_y; j++) {
            if (!(i >= in_st_x && i <= in_en_x && j >= in_st_y && j <= in_en_y))
                if (i > 0 && i < pic->cols && j > 0 && j < pic->rows) {
                    ext_val += pic->data[(pic->cols * j) + i];
                    val_cnt++;
                }
Thiago Santini's avatar
Thiago Santini committed
66 67
        }

68 69
    if (val_cnt > 0)
        ext_val = ext_val / val_cnt;
Thiago Santini's avatar
Thiago Santini committed
70 71 72
    else
        return false;

73
    val = ext_val - val;
Thiago Santini's avatar
Thiago Santini committed
74

75
    *erg = (int)val;
Thiago Santini's avatar
Thiago Santini committed
76

77 78 79 80
    if (val > 10)
        return true;
    else
        return false;
Thiago Santini's avatar
Thiago Santini committed
81 82
}

83 84
static int calc_inner_gray(Mat* pic, std::vector<Point> curve, RotatedRect ellipse)
{
Thiago Santini's avatar
Thiago Santini committed
85

86 87
    int gray_val = 0;
    int gray_cnt = 0;
Thiago Santini's avatar
Thiago Santini committed
88

89
    Mat checkmap = Mat::zeros(pic->size(), CV_8U);
Thiago Santini's avatar
Thiago Santini committed
90

91
    for (unsigned int i = 0; i < curve.size(); i++) {
Thiago Santini's avatar
Thiago Santini committed
92

93 94
        int vec_x = (int)round(curve[i].x - ellipse.center.x);
        int vec_y = (int)round(curve[i].y - ellipse.center.y);
Thiago Santini's avatar
Thiago Santini committed
95

96 97 98
        for (float p = 0.95f; p > 0.80f; p -= 0.01f) { //0.75;-0.05
            int p_x = (int)round(ellipse.center.x + float((float(vec_x) * p) + 0.5));
            int p_y = (int)round(ellipse.center.y + float((float(vec_y) * p) + 0.5));
Thiago Santini's avatar
Thiago Santini committed
99

100
            if (p_x > 0 && p_x < pic->cols && p_y > 0 && p_y < pic->rows) {
Thiago Santini's avatar
Thiago Santini committed
101

102 103 104
                if (checkmap.data[(pic->cols * p_y) + p_x] == 0) {
                    checkmap.data[(pic->cols * p_y) + p_x] = 1;
                    gray_val += (unsigned int)pic->data[(pic->cols * p_y) + p_x];
Thiago Santini's avatar
Thiago Santini committed
105 106 107 108 109 110
                    gray_cnt++;
                }
            }
        }
    }

111 112
    if (gray_cnt > 0)
        gray_val = gray_val / gray_cnt;
Thiago Santini's avatar
Thiago Santini committed
113
    else
114
        gray_val = 1000;
Thiago Santini's avatar
Thiago Santini committed
115 116 117 118

    return gray_val;
}

119 120
static std::vector<std::vector<Point>> get_curves(Mat* pic, Mat* edge, Mat* magni, int start_x, int end_x, int start_y, int end_y, double mean_dist, int inner_color_range)
{
Thiago Santini's avatar
Thiago Santini committed
121

122
    (void)magni;
Thiago Santini's avatar
Thiago Santini committed
123 124 125 126 127 128 129
    std::vector<std::vector<Point>> all_lines;

    std::vector<std::vector<Point>> all_curves;
    std::vector<Point> curve;

    std::vector<Point> all_means;

130 131 132 133 134 135 136 137
    if (start_x < 2)
        start_x = 2;
    if (start_y < 2)
        start_y = 2;
    if (end_x > pic->cols - 2)
        end_x = pic->cols - 2;
    if (end_y > pic->rows - 2)
        end_y = pic->rows - 2;
Thiago Santini's avatar
Thiago Santini committed
138

139
    int curve_idx = 0;
Thiago Santini's avatar
Thiago Santini committed
140 141 142
    Point mean_p;
    bool add_curve;
    int mean_inner_gray;
143
    int mean_inner_gray_last = 1000000;
Thiago Santini's avatar
Thiago Santini committed
144 145 146 147 148 149 150

    all_curves.clear();
    all_means.clear();
    all_lines.clear();

    bool check[IMG_SIZE][IMG_SIZE];

151 152 153
    for (int i = 0; i < IMG_SIZE; i++)
        for (int j = 0; j < IMG_SIZE; j++)
            check[i][j] = 0;
Thiago Santini's avatar
Thiago Santini committed
154 155

    //get all lines
156 157
    for (int i = start_x; i < end_x; i++)
        for (int j = start_y; j < end_y; j++) {
Thiago Santini's avatar
Thiago Santini committed
158

159 160
            if (edge->data[(edge->cols * (j)) + (i)] > 0 && !check[i][j]) {
                check[i][j] = 1;
Thiago Santini's avatar
Thiago Santini committed
161 162

                curve.clear();
163
                curve_idx = 0;
Thiago Santini's avatar
Thiago Santini committed
164

165 166 167
                curve.push_back(Point(i, j));
                mean_p.x = i;
                mean_p.y = j;
Thiago Santini's avatar
Thiago Santini committed
168 169
                curve_idx++;

170
                int akt_idx = 0;
Thiago Santini's avatar
Thiago Santini committed
171

172
                while (akt_idx < curve_idx) {
Thiago Santini's avatar
Thiago Santini committed
173

174 175 176
                    Point akt_pos = curve[akt_idx];
                    for (int k1 = -1; k1 < 2; k1++)
                        for (int k2 = -1; k2 < 2; k2++) {
Thiago Santini's avatar
Thiago Santini committed
177

178 179 180 181
                            if (akt_pos.x + k1 >= start_x && akt_pos.x + k1 < end_x && akt_pos.y + k2 >= start_y && akt_pos.y + k2 < end_y)
                                if (!check[akt_pos.x + k1][akt_pos.y + k2])
                                    if (edge->data[(edge->cols * (akt_pos.y + k2)) + (akt_pos.x + k1)] > 0) {
                                        check[akt_pos.x + k1][akt_pos.y + k2] = 1;
Thiago Santini's avatar
Thiago Santini committed
182

183 184 185 186 187 188
                                        mean_p.x += akt_pos.x + k1;
                                        mean_p.y += akt_pos.y + k2;
                                        curve.push_back(Point(akt_pos.x + k1, akt_pos.y + k2));
                                        curve_idx++;
                                    }
                        }
Thiago Santini's avatar
Thiago Santini committed
189 190 191
                    akt_idx++;
                }

192
                if (curve_idx > 10 && curve.size() > 10) {
Thiago Santini's avatar
Thiago Santini committed
193

194 195
                    mean_p.x = (int)floor((double(mean_p.x) / double(curve_idx)) + 0.5);
                    mean_p.y = (int)floor((double(mean_p.y) / double(curve_idx)) + 0.5);
Thiago Santini's avatar
Thiago Santini committed
196 197 198 199 200 201 202 203 204

                    all_means.push_back(mean_p);
                    all_lines.push_back(curve);
                }
            }
        }

    RotatedRect selected_ellipse;

205
    for (unsigned int iii = 0; iii < all_lines.size(); iii++) {
Thiago Santini's avatar
Thiago Santini committed
206

207 208
        curve = all_lines.at(iii);
        mean_p = all_means.at(iii);
Thiago Santini's avatar
Thiago Santini committed
209

210 211
        int ergebniss = 0;
        add_curve = true;
Thiago Santini's avatar
Thiago Santini committed
212 213 214

        RotatedRect ellipse;

215 216 217
        for (unsigned int i = 0; i < curve.size(); i++)
            if (abs(mean_p.x - curve[i].x) <= mean_dist && abs(mean_p.y - curve[i].y) <= mean_dist)
                add_curve = false;
Thiago Santini's avatar
Thiago Santini committed
218

219 220
        //is ellipse fit possible
        if (add_curve) {
Thiago Santini's avatar
Thiago Santini committed
221

222
            ellipse = fitEllipse(Mat(curve));
Thiago Santini's avatar
Thiago Santini committed
223

224
            if (ellipse.center.x < 0 || ellipse.center.y < 0 || ellipse.center.x > pic->cols || ellipse.center.y > pic->rows) {
Thiago Santini's avatar
Thiago Santini committed
225

226 227
                add_curve = false;
            }
Thiago Santini's avatar
Thiago Santini committed
228

229
            if (ellipse.size.height > 3 * ellipse.size.width || ellipse.size.width > 3 * ellipse.size.height) {
Thiago Santini's avatar
Thiago Santini committed
230

231 232
                add_curve = false;
            }
Thiago Santini's avatar
Thiago Santini committed
233

234 235 236 237
            if (add_curve) { // pupil area
                if (ellipse.size.width * ellipse.size.height < ElSe::minArea || ellipse.size.width * ellipse.size.height > ElSe::maxArea)
                    add_curve = false;
            }
Thiago Santini's avatar
Thiago Santini committed
238

239 240 241 242 243
            if (add_curve) {
                if (!is_good_ellipse_eval(&ellipse, pic, &ergebniss))
                    add_curve = false;
            }
        }
Thiago Santini's avatar
Thiago Santini committed
244

245
        if (add_curve) {
Thiago Santini's avatar
Thiago Santini committed
246

247 248
            if (inner_color_range >= 0) {
                mean_inner_gray = 0;
Thiago Santini's avatar
Thiago Santini committed
249

250 251
                mean_inner_gray = calc_inner_gray(pic, curve, ellipse);
                mean_inner_gray = (int)(mean_inner_gray * (1 + abs(ellipse.size.height - ellipse.size.width)));
Thiago Santini's avatar
Thiago Santini committed
252

253 254 255 256 257
                if (mean_inner_gray_last > mean_inner_gray) {
                    mean_inner_gray_last = mean_inner_gray;
                    all_curves.clear();
                    all_curves.push_back(curve);
                } else if (mean_inner_gray_last == mean_inner_gray) {
Thiago Santini's avatar
Thiago Santini committed
258

259 260 261 262 263
                    if (curve.size() > all_curves[0].size()) {
                        mean_inner_gray_last = mean_inner_gray;
                        all_curves.clear();
                        all_curves.push_back(curve);
                        selected_ellipse = ellipse;
Thiago Santini's avatar
Thiago Santini committed
264 265
                    }
                }
266
            }
Thiago Santini's avatar
Thiago Santini committed
267
        }
268
    }
Thiago Santini's avatar
Thiago Santini committed
269

270
    return all_curves;
Thiago Santini's avatar
Thiago Santini committed
271 272
}

273 274
static RotatedRect find_best_edge(Mat* pic, Mat* edge, Mat* magni, int start_x, int end_x, int start_y, int end_y, double mean_dist, int inner_color_range)
{
Thiago Santini's avatar
Thiago Santini committed
275 276

    RotatedRect ellipse;
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
    ellipse.center.x = 0;
    ellipse.center.y = 0;
    ellipse.angle = 0.0;
    ellipse.size.height = 0.0;
    ellipse.size.width = 0.0;

    std::vector<std::vector<Point>> all_curves = get_curves(pic, edge, magni, start_x, end_x, start_y, end_y, mean_dist, inner_color_range);

    if (all_curves.size() == 1) {
        ellipse = fitEllipse(Mat(all_curves[0]));

        if (ellipse.center.x < 0 || ellipse.center.y < 0 || ellipse.center.x > pic->cols || ellipse.center.y > pic->rows) {
            ellipse.center.x = 0;
            ellipse.center.y = 0;
            ellipse.angle = 0.0;
            ellipse.size.height = 0.0;
            ellipse.size.width = 0.0;
Thiago Santini's avatar
Thiago Santini committed
294 295
        }

296 297 298 299 300 301
    } else {
        ellipse.center.x = 0;
        ellipse.center.y = 0;
        ellipse.angle = 0.0;
        ellipse.size.height = 0.0;
        ellipse.size.width = 0.0;
Thiago Santini's avatar
Thiago Santini committed
302 303 304 305 306
    }

    return ellipse;
}

307 308
static void filter_edges(Mat* edge, int start_xx, int end_xx, int start_yy, int end_yy)
{
Thiago Santini's avatar
Thiago Santini committed
309

310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325
    int start_x = start_xx + 5;
    int end_x = end_xx - 5;
    int start_y = start_yy + 5;
    int end_y = end_yy - 5;

    if (start_x < 5)
        start_x = 5;
    if (end_x > edge->cols - 5)
        end_x = edge->cols - 5;
    if (start_y < 5)
        start_y = 5;
    if (end_y > edge->rows - 5)
        end_y = edge->rows - 5;

    for (int j = start_y; j < end_y; j++)
        for (int i = start_x; i < end_x; i++) {
Thiago Santini's avatar
Thiago Santini committed
326 327
            int box[9];

328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343
            box[4] = (int)edge->data[(edge->cols * (j)) + (i)];

            if (box[4]) {
                box[1] = (int)edge->data[(edge->cols * (j - 1)) + (i)];
                box[3] = (int)edge->data[(edge->cols * (j)) + (i - 1)];
                box[5] = (int)edge->data[(edge->cols * (j)) + (i + 1)];
                box[7] = (int)edge->data[(edge->cols * (j + 1)) + (i)];

                if ((box[5] && box[7]))
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if ((box[5] && box[1]))
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if ((box[3] && box[7]))
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if ((box[3] && box[1]))
                    edge->data[(edge->cols * (j)) + (i)] = 0;
Thiago Santini's avatar
Thiago Santini committed
344 345 346
            }
        }

347 348 349 350
    //too many neigbours
    for (int j = start_y; j < end_y; j++)
        for (int i = start_x; i < end_x; i++) {
            int neig = 0;
Thiago Santini's avatar
Thiago Santini committed
351

352 353
            for (int k1 = -1; k1 < 2; k1++)
                for (int k2 = -1; k2 < 2; k2++) {
Thiago Santini's avatar
Thiago Santini committed
354

355
                    if (edge->data[(edge->cols * (j + k1)) + (i + k2)] > 0)
Thiago Santini's avatar
Thiago Santini committed
356 357 358
                        neig++;
                }

359 360
            if (neig > 3)
                edge->data[(edge->cols * (j)) + (i)] = 0;
Thiago Santini's avatar
Thiago Santini committed
361 362
        }

363 364
    for (int j = start_y; j < end_y; j++)
        for (int i = start_x; i < end_x; i++) {
Thiago Santini's avatar
Thiago Santini committed
365 366
            int box[17];

367
            box[4] = (int)edge->data[(edge->cols * (j)) + (i)];
Thiago Santini's avatar
Thiago Santini committed
368

369 370 371 372
            if (box[4]) {
                box[0] = (int)edge->data[(edge->cols * (j - 1)) + (i - 1)];
                box[1] = (int)edge->data[(edge->cols * (j - 1)) + (i)];
                box[2] = (int)edge->data[(edge->cols * (j - 1)) + (i + 1)];
Thiago Santini's avatar
Thiago Santini committed
373

374 375
                box[3] = (int)edge->data[(edge->cols * (j)) + (i - 1)];
                box[5] = (int)edge->data[(edge->cols * (j)) + (i + 1)];
Thiago Santini's avatar
Thiago Santini committed
376

377 378 379
                box[6] = (int)edge->data[(edge->cols * (j + 1)) + (i - 1)];
                box[7] = (int)edge->data[(edge->cols * (j + 1)) + (i)];
                box[8] = (int)edge->data[(edge->cols * (j + 1)) + (i + 1)];
Thiago Santini's avatar
Thiago Santini committed
380 381

                //external
382 383
                box[9] = (int)edge->data[(edge->cols * (j)) + (i + 2)];
                box[10] = (int)edge->data[(edge->cols * (j + 2)) + (i)];
Thiago Santini's avatar
Thiago Santini committed
384

385 386 387
                box[11] = (int)edge->data[(edge->cols * (j)) + (i + 3)];
                box[12] = (int)edge->data[(edge->cols * (j - 1)) + (i + 2)];
                box[13] = (int)edge->data[(edge->cols * (j + 1)) + (i + 2)];
Thiago Santini's avatar
Thiago Santini committed
388

389 390 391
                box[14] = (int)edge->data[(edge->cols * (j + 3)) + (i)];
                box[15] = (int)edge->data[(edge->cols * (j + 2)) + (i - 1)];
                box[16] = (int)edge->data[(edge->cols * (j + 2)) + (i + 1)];
Thiago Santini's avatar
Thiago Santini committed
392

393 394 395 396
                if ((box[10] && !box[7]) && (box[8] || box[6])) {
                    edge->data[(edge->cols * (j + 1)) + (i - 1)] = 0;
                    edge->data[(edge->cols * (j + 1)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j + 1)) + (i)] = 255;
Thiago Santini's avatar
Thiago Santini committed
397 398
                }

399 400 401 402 403 404 405
                if ((box[14] && !box[7] && !box[10]) && ((box[8] || box[6]) && (box[16] || box[15]))) {
                    edge->data[(edge->cols * (j + 1)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j + 1)) + (i - 1)] = 0;
                    edge->data[(edge->cols * (j + 2)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j + 2)) + (i - 1)] = 0;
                    edge->data[(edge->cols * (j + 1)) + (i)] = 255;
                    edge->data[(edge->cols * (j + 2)) + (i)] = 255;
Thiago Santini's avatar
Thiago Santini committed
406 407
                }

408 409 410 411
                if ((box[9] && !box[5]) && (box[8] || box[2])) {
                    edge->data[(edge->cols * (j + 1)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j - 1)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j)) + (i + 1)] = 255;
Thiago Santini's avatar
Thiago Santini committed
412 413
                }

414 415 416 417 418 419 420
                if ((box[11] && !box[5] && !box[9]) && ((box[8] || box[2]) && (box[13] || box[12]))) {
                    edge->data[(edge->cols * (j + 1)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j - 1)) + (i + 1)] = 0;
                    edge->data[(edge->cols * (j + 1)) + (i + 2)] = 0;
                    edge->data[(edge->cols * (j - 1)) + (i + 2)] = 0;
                    edge->data[(edge->cols * (j)) + (i + 1)] = 255;
                    edge->data[(edge->cols * (j)) + (i + 2)] = 255;
Thiago Santini's avatar
Thiago Santini committed
421 422 423 424
                }
            }
        }

425 426
    for (int j = start_y; j < end_y; j++)
        for (int i = start_x; i < end_x; i++) {
Thiago Santini's avatar
Thiago Santini committed
427 428 429

            int box[33];

430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499
            box[4] = (int)edge->data[(edge->cols * (j)) + (i)];

            if (box[4]) {
                box[0] = (int)edge->data[(edge->cols * (j - 1)) + (i - 1)];
                box[1] = (int)edge->data[(edge->cols * (j - 1)) + (i)];
                box[2] = (int)edge->data[(edge->cols * (j - 1)) + (i + 1)];

                box[3] = (int)edge->data[(edge->cols * (j)) + (i - 1)];
                box[5] = (int)edge->data[(edge->cols * (j)) + (i + 1)];

                box[6] = (int)edge->data[(edge->cols * (j + 1)) + (i - 1)];
                box[7] = (int)edge->data[(edge->cols * (j + 1)) + (i)];
                box[8] = (int)edge->data[(edge->cols * (j + 1)) + (i + 1)];

                box[9] = (int)edge->data[(edge->cols * (j - 1)) + (i + 2)];
                box[10] = (int)edge->data[(edge->cols * (j - 1)) + (i - 2)];
                box[11] = (int)edge->data[(edge->cols * (j + 1)) + (i + 2)];
                box[12] = (int)edge->data[(edge->cols * (j + 1)) + (i - 2)];

                box[13] = (int)edge->data[(edge->cols * (j - 2)) + (i - 1)];
                box[14] = (int)edge->data[(edge->cols * (j - 2)) + (i + 1)];
                box[15] = (int)edge->data[(edge->cols * (j + 2)) + (i - 1)];
                box[16] = (int)edge->data[(edge->cols * (j + 2)) + (i + 1)];

                box[17] = (int)edge->data[(edge->cols * (j - 3)) + (i - 1)];
                box[18] = (int)edge->data[(edge->cols * (j - 3)) + (i + 1)];
                box[19] = (int)edge->data[(edge->cols * (j + 3)) + (i - 1)];
                box[20] = (int)edge->data[(edge->cols * (j + 3)) + (i + 1)];

                box[21] = (int)edge->data[(edge->cols * (j + 1)) + (i + 3)];
                box[22] = (int)edge->data[(edge->cols * (j + 1)) + (i - 3)];
                box[23] = (int)edge->data[(edge->cols * (j - 1)) + (i + 3)];
                box[24] = (int)edge->data[(edge->cols * (j - 1)) + (i - 3)];

                box[25] = (int)edge->data[(edge->cols * (j - 2)) + (i - 2)];
                box[26] = (int)edge->data[(edge->cols * (j + 2)) + (i + 2)];
                box[27] = (int)edge->data[(edge->cols * (j - 2)) + (i + 2)];
                box[28] = (int)edge->data[(edge->cols * (j + 2)) + (i - 2)];

                box[29] = (int)edge->data[(edge->cols * (j - 3)) + (i - 3)];
                box[30] = (int)edge->data[(edge->cols * (j + 3)) + (i + 3)];
                box[31] = (int)edge->data[(edge->cols * (j - 3)) + (i + 3)];
                box[32] = (int)edge->data[(edge->cols * (j + 3)) + (i - 3)];

                if (box[7] && box[2] && box[9])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[7] && box[0] && box[10])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[1] && box[8] && box[11])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[1] && box[6] && box[12])
                    edge->data[(edge->cols * (j)) + (i)] = 0;

                if (box[0] && box[13] && box[17] && box[8] && box[11] && box[21])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[2] && box[14] && box[18] && box[6] && box[12] && box[22])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[6] && box[15] && box[19] && box[2] && box[9] && box[23])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[8] && box[16] && box[20] && box[0] && box[10] && box[24])
                    edge->data[(edge->cols * (j)) + (i)] = 0;

                if (box[0] && box[25] && box[2] && box[27])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[0] && box[25] && box[6] && box[28])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[8] && box[26] && box[2] && box[27])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box[8] && box[26] && box[6] && box[28])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
Thiago Santini's avatar
Thiago Santini committed
500 501

                int box2[18];
502
                box2[1] = (int)edge->data[(edge->cols * (j)) + (i - 1)];
Thiago Santini's avatar
Thiago Santini committed
503

504 505
                box2[2] = (int)edge->data[(edge->cols * (j - 1)) + (i - 2)];
                box2[3] = (int)edge->data[(edge->cols * (j - 2)) + (i - 3)];
Thiago Santini's avatar
Thiago Santini committed
506

507 508
                box2[4] = (int)edge->data[(edge->cols * (j - 1)) + (i + 1)];
                box2[5] = (int)edge->data[(edge->cols * (j - 2)) + (i + 2)];
Thiago Santini's avatar
Thiago Santini committed
509

510 511
                box2[6] = (int)edge->data[(edge->cols * (j + 1)) + (i - 2)];
                box2[7] = (int)edge->data[(edge->cols * (j + 2)) + (i - 3)];
Thiago Santini's avatar
Thiago Santini committed
512

513 514
                box2[8] = (int)edge->data[(edge->cols * (j + 1)) + (i + 1)];
                box2[9] = (int)edge->data[(edge->cols * (j + 2)) + (i + 2)];
Thiago Santini's avatar
Thiago Santini committed
515

516
                box2[10] = (int)edge->data[(edge->cols * (j + 1)) + (i)];
Thiago Santini's avatar
Thiago Santini committed
517

518 519
                box2[15] = (int)edge->data[(edge->cols * (j - 1)) + (i - 1)];
                box2[16] = (int)edge->data[(edge->cols * (j - 2)) + (i - 2)];
Thiago Santini's avatar
Thiago Santini committed
520

521 522
                box2[11] = (int)edge->data[(edge->cols * (j + 2)) + (i + 1)];
                box2[12] = (int)edge->data[(edge->cols * (j + 3)) + (i + 2)];
Thiago Santini's avatar
Thiago Santini committed
523

524 525
                box2[13] = (int)edge->data[(edge->cols * (j + 2)) + (i - 1)];
                box2[14] = (int)edge->data[(edge->cols * (j + 3)) + (i - 2)];
Thiago Santini's avatar
Thiago Santini committed
526

527 528 529 530 531 532 533 534
                if (box2[1] && box2[2] && box2[3] && box2[4] && box2[5])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box2[1] && box2[6] && box2[7] && box2[8] && box2[9])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box2[10] && box2[11] && box2[12] && box2[4] && box2[5])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
                if (box2[10] && box2[13] && box2[14] && box2[15] && box2[16])
                    edge->data[(edge->cols * (j)) + (i)] = 0;
Thiago Santini's avatar
Thiago Santini committed
535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
            }
        }
}

//static float hypot(float a,float b)
//{
//    a=abs(a);
//    b=abs(b);
//    float t=a<b?a:b;
//    float x=a>b?a:b;
//
//    if (x==0)
//        return 0;
//
//    t = t/x;
//    return x*sqrt(1+(t*t));
//}

#define MAX_LINE 10000

555 556
static void matlab_bwselect(Mat* strong, Mat* weak, Mat* check)
{
Thiago Santini's avatar
Thiago Santini committed
557

558 559
    int pic_x = strong->cols;
    int pic_y = strong->rows;
Thiago Santini's avatar
Thiago Santini committed
560 561

    int lines[MAX_LINE];
562
    int lines_idx = 0;
Thiago Santini's avatar
Thiago Santini committed
563

564
    int idx = 0;
Thiago Santini's avatar
Thiago Santini committed
565

566
    for (int i = 1; i < pic_y - 1; i++) {
Thiago Santini's avatar
Thiago Santini committed
567

568
        for (int j = 1; j < pic_x - 1; j++) {
Thiago Santini's avatar
Thiago Santini committed
569

570
            if (strong->data[idx + j] != 0 && check->data[idx + j] == 0) {
Thiago Santini's avatar
Thiago Santini committed
571

572 573 574
                check->data[idx + j] = 255;
                lines_idx = 1;
                lines[0] = idx + j;
Thiago Santini's avatar
Thiago Santini committed
575

576
                int akt_idx = 0;
Thiago Santini's avatar
Thiago Santini committed
577

578
                while (akt_idx < lines_idx && lines_idx < MAX_LINE - 1) {
Thiago Santini's avatar
Thiago Santini committed
579

580
                    int akt_pos = lines[akt_idx];
Thiago Santini's avatar
Thiago Santini committed
581

582 583 584
                    if (akt_pos - pic_x - 1 >= 0 && akt_pos + pic_x + 1 < pic_x * pic_y) {
                        for (int k1 = -1; k1 < 2; k1++)
                            for (int k2 = -1; k2 < 2; k2++) {
Thiago Santini's avatar
Thiago Santini committed
585

586
                                if (check->data[(akt_pos + (k1 * pic_x)) + k2] == 0 && weak->data[(akt_pos + (k1 * pic_x)) + k2] != 0) {
Thiago Santini's avatar
Thiago Santini committed
587

588
                                    check->data[(akt_pos + (k1 * pic_x)) + k2] = 255;
Thiago Santini's avatar
Thiago Santini committed
589

590 591 592
                                    lines_idx++;
                                    lines[lines_idx - 1] = (akt_pos + (k1 * pic_x)) + k2;
                                }
Thiago Santini's avatar
Thiago Santini committed
593 594 595 596 597 598 599
                            }
                    }
                    akt_idx++;
                }
            }
        }

600
        idx += pic_x;
Thiago Santini's avatar
Thiago Santini committed
601 602 603
    }
}

604 605 606
static Mat canny_impl(Mat* pic, Mat* magni)
{
    int k_sz = 16;
Thiago Santini's avatar
Thiago Santini committed
607

608 609 610 611 612 613 614
    float gau[16] = { 0.000000220358050f, 0.000007297256405f, 0.000146569312970f, 0.001785579770079f,
        0.013193749090229f, 0.059130281094460f, 0.160732768610747f, 0.265003534507060f, 0.265003534507060f,
        0.160732768610747f, 0.059130281094460f, 0.013193749090229f, 0.001785579770079f, 0.000146569312970f,
        0.000007297256405f, 0.000000220358050f };
    float deriv_gau[16] = { -0.000026704586264f, -0.000276122963398f, -0.003355163265098f, -0.024616683775044f, -0.108194751875585f,
        -0.278368310241814f, -0.388430056419619f, -0.196732206873178f, 0.196732206873178f, 0.388430056419619f,
        0.278368310241814f, 0.108194751875585f, 0.024616683775044f, 0.003355163265098f, 0.000276122963398f, 0.000026704586264f };
Thiago Santini's avatar
Thiago Santini committed
615

616
    Point anchor = Point(-1, -1);
Thiago Santini's avatar
Thiago Santini committed
617 618 619 620 621
    float delta = 0;
    int ddepth = -1;

    pic->convertTo(*pic, CV_32FC1);

622 623
    Mat gau_x = Mat(1, k_sz, CV_32FC1, &gau);
    Mat deriv_gau_x = Mat(1, k_sz, CV_32FC1, &deriv_gau);
Thiago Santini's avatar
Thiago Santini committed
624 625 626 627

    Mat res_x;
    Mat res_y;

628 629 630 631
    transpose(*pic, *pic);
    filter2D(*pic, res_x, ddepth, gau_x, anchor, delta, BORDER_REPLICATE);
    transpose(*pic, *pic);
    transpose(res_x, res_x);
Thiago Santini's avatar
Thiago Santini committed
632

633
    filter2D(res_x, res_x, ddepth, deriv_gau_x, anchor, delta, BORDER_REPLICATE);
Thiago Santini's avatar
Thiago Santini committed
634

635
    filter2D(*pic, res_y, ddepth, gau_x, anchor, delta, BORDER_REPLICATE);
Thiago Santini's avatar
Thiago Santini committed
636

637 638 639
    transpose(res_y, res_y);
    filter2D(res_y, res_y, ddepth, deriv_gau_x, anchor, delta, BORDER_REPLICATE);
    transpose(res_y, res_y);
Thiago Santini's avatar
Thiago Santini committed
640

641
    *magni = Mat::zeros(pic->rows, pic->cols, CV_32FC1);
Thiago Santini's avatar
Thiago Santini committed
642

643 644 645 646 647
    float *p_res, *p_x, *p_y;
    for (int i = 0; i < magni->rows; i++) {
        p_res = magni->ptr<float>(i);
        p_x = res_x.ptr<float>(i);
        p_y = res_y.ptr<float>(i);
Thiago Santini's avatar
Thiago Santini committed
648

649
        for (int j = 0; j < magni->cols; j++) {
Thiago Santini's avatar
Thiago Santini committed
650 651 652 653 654 655
            //res.at<float>(j, i)= sqrt( (res_x.at<float>(j, i)*res_x.at<float>(j, i)) + (res_y.at<float>(j, i)*res_y.at<float>(j, i)) );
            //res.at<float>(j, i)=robust_pytagoras_after_MOLAR_MORRIS(res_x.at<float>(j, i), res_y.at<float>(j, i));
            //res.at<float>(j, i)=hypot(res_x.at<float>(j, i), res_y.at<float>(j, i));

            //p_res[j]=__ieee754_hypot(p_x[j], p_y[j]);

656
            p_res[j] = hypot(p_x[j], p_y[j]);
Thiago Santini's avatar
Thiago Santini committed
657 658 659 660 661
        }
    }

    //th selection

662 663
    int PercentOfPixelsNotEdges = (int)round(0.7 * magni->cols * magni->rows);
    float ThresholdRatio = 0.4f;
Thiago Santini's avatar
Thiago Santini committed
664

665 666
    float high_th = 0;
    float low_th = 0;
Thiago Santini's avatar
Thiago Santini committed
667

668
    int h_sz = 64;
Thiago Santini's avatar
Thiago Santini committed
669
    int hist[64];
670 671
    for (int i = 0; i < h_sz; i++)
        hist[i] = 0;
Thiago Santini's avatar
Thiago Santini committed
672 673 674

    normalize(*magni, *magni, 0, 1, NORM_MINMAX, CV_32FC1);

675
    Mat res_idx = Mat::zeros(pic->rows, pic->cols, CV_8U);
Thiago Santini's avatar
Thiago Santini committed
676 677
    normalize(*magni, res_idx, 0, 63, NORM_MINMAX, CV_32S);

678 679 680 681
    int* p_res_idx = 0;
    for (int i = 0; i < magni->rows; i++) {
        p_res_idx = res_idx.ptr<int>(i);
        for (int j = 0; j < magni->cols; j++) {
Thiago Santini's avatar
Thiago Santini committed
682
            hist[p_res_idx[j]]++;
683 684
        }
    }
Thiago Santini's avatar
Thiago Santini committed
685

686 687 688 689 690
    int sum = 0;
    for (int i = 0; i < h_sz; i++) {
        sum += hist[i];
        if (sum > PercentOfPixelsNotEdges) {
            high_th = float(i + 1) / float(h_sz);
Thiago Santini's avatar
Thiago Santini committed
691 692 693
            break;
        }
    }
694
    low_th = ThresholdRatio * high_th;
Thiago Santini's avatar
Thiago Santini committed
695 696

    //non maximum supression + interpolation
697 698
    Mat non_ms = Mat::zeros(pic->rows, pic->cols, CV_8U);
    Mat non_ms_hth = Mat::zeros(pic->rows, pic->cols, CV_8U);
Thiago Santini's avatar
Thiago Santini committed
699

700
    float ix, iy, grad1, grad2, d;
Thiago Santini's avatar
Thiago Santini committed
701

702 703 704 705 706
    char *p_non_ms, *p_non_ms_hth;
    float *p_res_t, *p_res_b;
    for (int i = 1; i < magni->rows - 1; i++) {
        p_non_ms = non_ms.ptr<char>(i);
        p_non_ms_hth = non_ms_hth.ptr<char>(i);
Thiago Santini's avatar
Thiago Santini committed
707

708 709 710
        p_res = magni->ptr<float>(i);
        p_res_t = magni->ptr<float>(i - 1);
        p_res_b = magni->ptr<float>(i + 1);
Thiago Santini's avatar
Thiago Santini committed
711

712 713
        p_x = res_x.ptr<float>(i);
        p_y = res_y.ptr<float>(i);
Thiago Santini's avatar
Thiago Santini committed
714

715
        for (int j = 1; j < magni->cols - 1; j++) {
Thiago Santini's avatar
Thiago Santini committed
716

717 718
            iy = p_y[j];
            ix = p_x[j];
Thiago Santini's avatar
Thiago Santini committed
719

720
            if ((iy <= 0 && ix > -iy) || (iy >= 0 && ix < -iy)) {
Thiago Santini's avatar
Thiago Santini committed
721

722 723 724
                d = abs(iy / ix);
                grad1 = (p_res[j + 1] * (1 - d)) + (p_res_t[j + 1] * d);
                grad2 = (p_res[j - 1] * (1 - d)) + (p_res_b[j - 1] * d);
Thiago Santini's avatar
Thiago Santini committed
725

726 727
                if (p_res[j] >= grad1 && p_res[j] >= grad2) {
                    p_non_ms[j] = (char)255;
Thiago Santini's avatar
Thiago Santini committed
728

729 730
                    if (p_res[j] > high_th)
                        p_non_ms_hth[j] = (char)255;
Thiago Santini's avatar
Thiago Santini committed
731
                }
732
            }
Thiago Santini's avatar
Thiago Santini committed
733

734 735 736 737
            if ((ix > 0 && -iy >= ix) || (ix < 0 && -iy <= ix)) {
                d = abs(ix / iy);
                grad1 = (p_res_t[j] * (1 - d)) + (p_res_t[j + 1] * d);
                grad2 = (p_res_b[j] * (1 - d)) + (p_res_b[j - 1] * d);
Thiago Santini's avatar
Thiago Santini committed
738

739 740 741 742
                if (p_res[j] >= grad1 && p_res[j] >= grad2) {
                    p_non_ms[j] = (char)255;
                    if (p_res[j] > high_th)
                        p_non_ms_hth[j] = (char)255;
Thiago Santini's avatar
Thiago Santini committed
743
                }
744
            }
Thiago Santini's avatar
Thiago Santini committed
745

746 747 748 749
            if ((ix <= 0 && ix > iy) || (ix >= 0 && ix < iy)) {
                d = abs(ix / iy);
                grad1 = (p_res_t[j] * (1 - d)) + (p_res_t[j - 1] * d);
                grad2 = (p_res_b[j] * (1 - d)) + (p_res_b[j + 1] * d);
Thiago Santini's avatar
Thiago Santini committed
750

751 752 753 754
                if (p_res[j] >= grad1 && p_res[j] >= grad2) {
                    p_non_ms[j] = (char)255;
                    if (p_res[j] > high_th)
                        p_non_ms_hth[j] = (char)255;
Thiago Santini's avatar
Thiago Santini committed
755
                }
756
            }
Thiago Santini's avatar
Thiago Santini committed
757

758 759 760 761
            if ((iy < 0 && ix <= iy) || (iy > 0 && ix >= iy)) {
                d = abs(iy / ix);
                grad1 = (p_res[j - 1] * (1 - d)) + (p_res_t[j - 1] * d);
                grad2 = (p_res[j + 1] * (1 - d)) + (p_res_b[j + 1] * d);
Thiago Santini's avatar
Thiago Santini committed
762

763 764 765 766
                if (p_res[j] >= grad1 && p_res[j] >= grad2) {
                    p_non_ms[j] = (char)255;
                    if (p_res[j] > high_th)
                        p_non_ms_hth[j] = (char)255;
Thiago Santini's avatar
Thiago Santini committed
767
                }
768 769 770
            }
        }
    }
Thiago Santini's avatar
Thiago Santini committed
771 772

    ////bw select
773 774
    Mat res_lin = Mat::zeros(pic->rows, pic->cols, CV_8U);
    matlab_bwselect(&non_ms_hth, &non_ms, &res_lin);
Thiago Santini's avatar
Thiago Santini committed
775 776 777 778 779 780

    pic->convertTo(*pic, CV_8U);

    return res_lin;
}

781 782
static void mum(Mat* pic, Mat* result, int fak)
{
Thiago Santini's avatar
Thiago Santini committed
783

784 785 786
    int fak_ges = fak + 1;
    int sz_x = pic->cols / fak_ges;
    int sz_y = pic->rows / fak_ges;
Thiago Santini's avatar
Thiago Santini committed
787

788
    *result = Mat::zeros(sz_y, sz_x, CV_8U);
Thiago Santini's avatar
Thiago Santini committed
789 790

    int hist[256];
791 792 793
    int mean = 0;
    int cnt = 0;
    int mean_2 = 0;
Thiago Santini's avatar
Thiago Santini committed
794

795 796
    int idx = 0;
    int idy = 0;
Thiago Santini's avatar
Thiago Santini committed
797

798 799
    for (int i = 0; i < sz_y; i++) {
        idy += fak_ges;
Thiago Santini's avatar
Thiago Santini committed
800

801 802
        for (int j = 0; j < sz_x; j++) {
            idx += fak_ges;
Thiago Santini's avatar
Thiago Santini committed
803

804 805
            for (int k = 0; k < 256; k++)
                hist[k] = 0;
Thiago Santini's avatar
Thiago Santini committed
806

807 808
            mean = 0;
            cnt = 0;
Thiago Santini's avatar
Thiago Santini committed
809

810 811
            for (int ii = -fak; ii <= fak; ii++)
                for (int jj = -fak; jj <= fak; jj++) {
Thiago Santini's avatar
Thiago Santini committed
812

813 814 815
                    if (idy + ii > 0 && idy + ii < pic->rows && idx + jj > 0 && idx + jj < pic->cols) {
                        if ((unsigned int)pic->data[(pic->cols * (idy + ii)) + (idx + jj)] > 255)
                            pic->data[(pic->cols * (idy + ii)) + (idx + jj)] = 255;
Thiago Santini's avatar
Thiago Santini committed
816

817
                        hist[pic->data[(pic->cols * (idy + ii)) + (idx + jj)]]++;
Thiago Santini's avatar
Thiago Santini committed
818
                        cnt++;
819
                        mean += pic->data[(pic->cols * (idy + ii)) + (idx + jj)];
Thiago Santini's avatar
Thiago Santini committed
820 821 822
                    }
                }

823
            mean = mean / cnt;
Thiago Santini's avatar
Thiago Santini committed
824

825 826 827 828 829
            mean_2 = 0;
            cnt = 0;
            for (int ii = 0; ii <= mean; ii++) {
                mean_2 += ii * hist[ii];
                cnt += hist[ii];
Thiago Santini's avatar
Thiago Santini committed
830 831
            }

832 833
            if (cnt == 0)
                mean_2 = mean;
Thiago Santini's avatar
Thiago Santini committed
834
            else
835
                mean_2 = mean_2 / cnt;
Thiago Santini's avatar
Thiago Santini committed
836

837
            result->data[(sz_x * (i)) + (j)] = mean_2;
Thiago Santini's avatar
Thiago Santini committed
838 839
        }

840
        idx = 0;
Thiago Santini's avatar
Thiago Santini committed
841 842 843
    }
}

844 845
static void gen_blob_neu(int rad, Mat* all_mat, Mat* all_mat_neg)
{
Thiago Santini's avatar
Thiago Santini committed
846

847 848 849 850
    int len = 1 + (4 * rad);
    int c0 = rad * 2;
    float negis = 0;
    float posis = 0;
Thiago Santini's avatar
Thiago Santini committed
851

852 853
    *all_mat = Mat::zeros(len, len, CV_32FC1);
    *all_mat_neg = Mat::zeros(len, len, CV_32FC1);
Thiago Santini's avatar
Thiago Santini committed
854

855 856 857
    float *p, *p_neg;
    for (int i = -rad * 2; i <= rad * 2; i++) { //height
        p = all_mat->ptr<float>(c0 + i);
Thiago Santini's avatar
Thiago Santini committed
858

859
        for (int j = -rad * 2; j <= rad * 2; j++) {
Thiago Santini's avatar
Thiago Santini committed
860

861 862 863
            if (i < -rad || i > rad) { //pos
                p[c0 + j] = 1;
                posis++;
Thiago Santini's avatar
Thiago Santini committed
864

865
            } else { //neg
Thiago Santini's avatar
Thiago Santini committed
866

867
                int sz_w = (int)sqrt(float(rad * rad) - float(i * i));
Thiago Santini's avatar
Thiago Santini committed
868

869 870 871 872 873 874
                if (abs(j) <= sz_w) {
                    p[c0 + j] = -1;
                    negis++;
                } else {
                    p[c0 + j] = 1;
                    posis++;
Thiago Santini's avatar
Thiago Santini committed
875 876 877
                }
            }
        }
878
    }
Thiago Santini's avatar
Thiago Santini committed
879

880 881 882
    for (int i = 0; i < len; i++) {
        p = all_mat->ptr<float>(i);
        p_neg = all_mat_neg->ptr<float>(i);
Thiago Santini's avatar