findzc.c

Go to the documentation of this file.

#include <stdio.h>
#include <math.h>


#define TINY    1.0e-3


int
G_math_findzc(double conv[], int size, double zc[], double thresh,
              int num_orients)
{
    int i, j, p;

    /* go through entire conv image - but skip border rows and cols */
    for (i = 1; i < size - 1; i++) {
        for (p = i * size + 1, j = 1; j < size - 1; j++, p++) {
            int nbr[4];
            int ni;

            /* examine the 4-neighbors of position p */
            nbr[0] = p - 1;     /* left */
            nbr[1] = p + 1;     /* right */
            nbr[2] = p - size;  /* up */
            nbr[3] = p + size;  /* down */

            zc[p] = 0;

            for (ni = 0; ni < 4; ni++) {
                /* condition for a zc: sign is different than a neighbor
                 * and the absolute value is less than that neighbor.
                 * Also, threshold magnitudes to eliminate noise
                 */
                if ((((conv[p] > 0) && (conv[nbr[ni]] < 0)) ||
                     ((conv[p] < 0) && (conv[nbr[ni]] > 0))) &&
                    (fabs(conv[p]) < fabs(conv[nbr[ni]])) &&
                    (fabs(conv[p] - conv[nbr[ni]]) > thresh)) {
                    double ang;
                    int dir;

                    /* found a zc here, get angle of gradient */
                    if (fabs(conv[nbr[1]] - conv[nbr[0]]) < TINY) {
                        ang = M_PI_2;

                        if (conv[nbr[2]] - conv[nbr[3]] < 0)
                            ang = -ang;
                    }
                    else
                        ang = atan2(conv[nbr[2]] - conv[nbr[3]],
                                    conv[nbr[1]] - conv[nbr[0]]);

                    /* scale -PI..PI to 0..num_orients - 1 */
                    dir =
                        num_orients * ((ang + M_PI) / (M_PI * 2.0)) + 0.4999;

                    /* shift scale so that 0 (not 8) is straight down */
                    dir = (3 * num_orients / 4 + dir) % num_orients;

                    /* add to differentiate between no zc and an orientation */
                    zc[p] = 1 + dir;
                    break;      /* quit looking at neighbors */
                }
            }                   /* for ni */
        }                       /* for p */
    }

    return 0;
}