view release on metacpan or  search on metacpan

RectanglesContainingDot_XS.xs  view on Meta::CPAN

    current = algo->current;
    if (!current || (current->top >= RECTANGLES_CHUNK_SIZE))
        current = allocate_rectangles_chunk(aTHX_ algo);

    rect = &(current->rects[current->top++]);

    rect->x0 = x0;
    rect->y0 = y0;
    rect->x1 = x1;
    rect->y1 = y1;

    rect->name = newSVsv(name);
}

int
double_cmp(pTHX_ double *a, double *b) {
    double fa = *a;
    double fb = *b;
    /* printf("cmp(%f, %f) => %d\n", fa, fb, (fa < fb) ? -1 : ((fa > fb) ? 1 : 0));  */
    return (fa < fb) ? -1 : ((fa > fb) ? 1 : 0);
}

void
sort_inplace(pTHX_ double **v, int size) {
    sortsv((SV**)v, size, (SVCOMPARE_t)&double_cmp);
}

struct division *
init_division(pTHX_ struct algorithm *algo) {
    struct rectangles_chunk *chunk;
    struct division *div;
    struct rectangle **rect;
    int size, i;

    if (algo->div)
        return algo->div;

    /* printf("."); fflush(stdout); */
    
    size = 0;
    for(chunk = algo->chunks; chunk; chunk = chunk->next)
        size += chunk->top;

    div = allocate_division(aTHX_ size);

    rect = div->rects;
    for (chunk = algo->chunks; chunk; chunk = chunk->next) {
        int i;
        int top = chunk->top;
        for (i=0; i<top; i++, rect++) {
            /* printf("."); fflush(stdout); */
            *rect = &(chunk->rects[i]);
        }
    }
    
    return algo->div = div;
}

double
find_best_cut(pTHX_ struct rectangle **rects, int size, int dir,
              double *bestv, int *sizel, int *sizer) {
    double **v0, **v1, **vc0, **vc1;
    double v, med, best;
    int op, cl;
    int i;

    my_assert(bestv);
    my_assert(sizel);
    my_assert(sizer);
    
    Newy(v0, size + 1, double *);
    Newy(v1, size + 1, double *);
    
    v0[size] = v1[size] = NULL;
    
    vc0 = v0; vc1 = v1;
    
    for (i=0; i<size; i++) {
        if (dir == 'x') {
            v0[i] = &(rects[i]->x0);
            v1[i] = &(rects[i]->x1);
        }
        else {
            v0[i] = &(rects[i]->y0);
            v1[i] = &(rects[i]->y1);
        }
        /* printf( "%d: v0=%g, v1=%g\n", i, *(v0[i]), *(v1[i])); fflush(stdout); */
    }
    v0[size] = v1[size] = NULL;

    sort_inplace(aTHX_ v0, size);
    sort_inplace(aTHX_ v1, size);
    
    op = cl = 0;
    med = 0.5 * size;
    best = (double)size * (double)size;

    my_assert(best >= 0);
             
    
    while (*v0 && *v1) {
        double v, good;
        double l, r;
        
        v =  (**v0 <= **v1) ? **v0 : **v1;

        while (*v0 && v == **v0) {
            op++;
            v0++;
        }
        while (*v1 && v == **v1) {
            cl++;
            v1++;
        }

        my_assert(op > 0 && op <= size);
        my_assert(cl >= 0 && cl <= size);
        
        l = op - med;
        r = size - cl - med;
        good = (double)l * (double)l + (double)r * (double)r;

        my_assert(good >= 0);
        
        if (good < best) {
            best = good;
            *bestv = v;
            *sizel = op;
            *sizer = size - cl;
        }
    }

    Safefry(vc0);
    Safefry(vc1);
    
    return best;
}

void
part_division(pTHX_ struct rectangle **rects, int size,
              double cut, int dir,
              struct division **left, int left_size,
              struct division **right, int right_size) {

    int i;
    struct rectangle **rectsl, **rectsr;

    my_assert(left);
    my_assert(right);
    my_assert(left_size);
    my_assert(right_size);
    my_assert(right_size < size);
    my_assert(left_size < size);
    
    *left = allocate_division(aTHX_ left_size);
    *right = allocate_division(aTHX_ right_size);

    // fprintf(stderr, "%d => %d, %d\n", size, left_size, right_size);
    
    rectsl = (*left)->rects;
    rectsr = (*right)->rects;

    for (i = 0; i < size; i++) {
        struct rectangle *rect = rects[i];
        if (dir == 'x') {
            if (cut >= rect->x0) *(rectsl++) = rect;
            if (cut < rect->x1) *(rectsr++) = rect;
        }
        else {
            if (cut >= rect->y0) *(rectsl++) = rect;
            if (cut < rect->y1) *(rectsr++) = rect;
        }
    }
    my_assert(rectsl == (*left)->rects + (*left)->size);
    my_assert(rectsr == (*right)->rects + (*right)->size);
}

int
subdivide_division(pTHX_ struct division *div) {
    int size;

    my_assert(div);
    
    size = div->size;
    if (size > MIN_DIVISION) {
        struct rectangle **rects = div->rects;
        double bestreq = 0.24 * size * size;
        double bestx, bestxx, besty, bestyy;
        int sizelx, sizerx, sizely, sizery;
        
        bestx = find_best_cut(aTHX_ rects, size, 'x', &bestxx, &sizelx, &sizerx);

        if (bestx > 0)
            besty = find_best_cut(aTHX_ rects, size, 'y', &bestyy, &sizely, &sizery);
        else
            besty = 1;

        if (bestx < besty) {
            if (bestx < bestreq) {
                // fprintf(stderr, "bestx: %f, bestreq: %f\n", bestx, bestreq);
                part_division(aTHX_ rects, size, bestxx, 'x', &(div->left), sizelx, &(div->right), sizerx);
                div->cut = bestxx;
                Safefry(div->rects);
                div->rects = NULL;
                return div->dir = 'x';
            }
        }
        else {
            if (besty < bestreq) {
                // fprintf(stderr, "besty: %f, bestreq: %f\n", besty, bestreq);
                part_division(aTHX_ rects, size, bestyy, 'y', &(div->left), sizely, &(div->right), sizery);
                div->cut = bestyy;
                Safefry(div->rects);
                div->rects = NULL;
                return div->dir = 'y';
            }
        }
    }
    return div->dir = 'n';
}

struct division *
division_containing_dot(pTHX_ struct division *div, double x, double y) {
    while(1) {
        int dir = div->dir;
        if (!dir)
            dir = subdivide_division(aTHX_ div);

        switch(dir) {
        case 'x':
            div = (x <= div->cut) ? div->left : div->right;
            break;
        case 'y':
            div = (y <= div->cut) ? div->left : div->right;
            break;
        default:
            my_assert(div->rects);
            return div;
        }
    }
}



MODULE = Algorithm::RectanglesContainingDot_XS		PACKAGE = Algorithm::RectanglesContainingDot_XS		

void
add_rectangle(self, name, x0, y0, x1, y1)
    struct algorithm *self
    SV *name
    double x0
    double y0
    double x1
    double y1
C_ARGS:
    aTHX_ self, name, x0, y0, x1, y1
    
void
rectangles_containing_dot(self, x, y)
    struct algorithm *self
    double x
    double y
PREINIT:
    struct division *div;
    int n = 0;
PPCODE:
    div = init_division(aTHX_ self);
    if (div) {
        struct rectangle **rects;
        int i, size;
        
        div = division_containing_dot(aTHX_ div, x, y);