rt_raster_intersects_algorithm()

static int rt_raster_intersects_algorithm ( rt_raster  rast1, rt_raster  rast2, rt_band  band1, rt_band  band2, int  hasnodata1, int  hasnodata2, double  nodata1, double  nodata2  )

static

Definition at line 637 of file rt_spatial_relationship.c.

References ES_NONE, FALSE, FLT_EQ, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_pixel(), rt_raster_cell_to_geopoint(), rt_raster_geopoint_to_cell(), rt_raster_get_height(), rt_raster_get_width(), rt_raster_get_x_scale(), and rt_raster_get_y_scale().

Referenced by rt_raster_intersects().

  {
        int i;
        int byHeight = 1;
        uint32_t dimValue;

        uint32_t row;
        uint32_t rowoffset;
        uint32_t col;
        uint32_t coloffset;

        enum line_points {X1, Y1, X2, Y2};
        enum point {pX, pY};
        double line1[4] = {0.};
        double line2[4] = {0.};
        double P[2] = {0.};
        double Qw[2] = {0.};
        double Qr[2] = {0.};
        double gt1[6] = {0.};
        double gt2[6] = {0.};
        double igt1[6] = {0};
        double igt2[6] = {0};
        double d;
        double val1;
        int noval1;
        int isnodata1;
        double val2;
        int noval2;
        int isnodata2;
        uint32_t adjacent[8] = {0};

        double xscale;
        double yscale;

        uint16_t width1;
        uint16_t height1;
        uint16_t width2;
        uint16_t height2;

        width1 = rt_raster_get_width(rast1);
        height1 = rt_raster_get_height(rast1);
        width2 = rt_raster_get_width(rast2);
        height2 = rt_raster_get_height(rast2);

        /* sampling scale */
        xscale = fmin(rt_raster_get_x_scale(rast1), rt_raster_get_x_scale(rast2)) / 10.;
        yscale = fmin(rt_raster_get_y_scale(rast1), rt_raster_get_y_scale(rast2)) / 10.;

        /* see if skew made rast2's rows are parallel to rast1's cols */
        rt_raster_cell_to_geopoint(
                rast1,
                0, 0,
                &(line1[X1]), &(line1[Y1]),
                gt1
        );

        rt_raster_cell_to_geopoint(
                rast1,
                0, height1,
                &(line1[X2]), &(line1[Y2]),
                gt1
        );

        rt_raster_cell_to_geopoint(
                rast2,
                0, 0,
                &(line2[X1]), &(line2[Y1]),
                gt2
        );

        rt_raster_cell_to_geopoint(
                rast2,
                width2, 0,
                &(line2[X2]), &(line2[Y2]),
                gt2
        );

        /* Parallel lines */
        if (FLT_EQ(((line1[X2] - line1[X1]) * (line2[Y2] - line2[Y1])),
                   ((line2[X2] - line2[X1]) * (line1[Y2] - line1[Y1]))))
                byHeight = 0;

        if (byHeight)
                dimValue = height2;
        else
                dimValue = width2;
        RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);

        /* 3 x 3 search */
        for (coloffset = 0; coloffset < 3; coloffset++) {
                for (rowoffset = 0; rowoffset < 3; rowoffset++) {
                        /* smaller raster */
                        for (col = coloffset; col <= width1; col += 3) {

                                rt_raster_cell_to_geopoint(
                                        rast1,
                                        col, 0,
                                        &(line1[X1]), &(line1[Y1]),
                                        gt1
                                );

                                rt_raster_cell_to_geopoint(
                                        rast1,
                                        col, height1,
                                        &(line1[X2]), &(line1[Y2]),
                                        gt1
                                );

                                /* larger raster */
                                for (row = rowoffset; row <= dimValue; row += 3) {

                                        if (byHeight) {
                                                rt_raster_cell_to_geopoint(
                                                        rast2,
                                                        0, row,
                                                        &(line2[X1]), &(line2[Y1]),
                                                        gt2
                                                );

                                                rt_raster_cell_to_geopoint(
                                                        rast2,
                                                        width2, row,
                                                        &(line2[X2]), &(line2[Y2]),
                                                        gt2
                                                );
                                        }
                                        else {
                                                rt_raster_cell_to_geopoint(
                                                        rast2,
                                                        row, 0,
                                                        &(line2[X1]), &(line2[Y1]),
                                                        gt2
                                                );

                                                rt_raster_cell_to_geopoint(
                                                        rast2,
                                                        row, height2,
                                                        &(line2[X2]), &(line2[Y2]),
                                                        gt2
                                                );
                                        }

                                        RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
                                        RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
                                                line1[X1], line1[Y1], line1[X2], line1[Y2]);
                                        RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
                                                line2[X1], line2[Y1], line2[X2], line2[Y2]);

                                        /* intersection */
                                        /* http://en.wikipedia.org/wiki/Line-line_intersection */
                                        d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
                                        /* no intersection */
                                        if (FLT_EQ(d, 0.)) {
                                                continue;
                                        }

                                        P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
                                                ((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
                                        P[pX] = P[pX] / d;

                                        P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
                                                ((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
                                        P[pY] = P[pY] / d;

                                        RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);

                                        /* intersection within bounds */
                                        if ((
                                                        (FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
                                                                (P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
                                                ) && (
                                                        (FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
                                                                (P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
                                                ) && (
                                                        (FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
                                                                (P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
                                                ) && (
                                                        (FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
                                                                (P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
                                        )) {
                                                RASTER_DEBUG(4, "within bounds");

                                                for (i = 0; i < 8; i++) adjacent[i] = 0;

                                                /* test points around intersection */
                                                for (i = 0; i < 8; i++) {
                                                        switch (i) {
                                                                case 7:
                                                                        Qw[pX] = P[pX] - xscale;
                                                                        Qw[pY] = P[pY] + yscale;
                                                                        break;
                                                                /* 270 degrees = 09:00 */
                                                                case 6:
                                                                        Qw[pX] = P[pX] - xscale;
                                                                        Qw[pY] = P[pY];
                                                                        break;
                                                                case 5:
                                                                        Qw[pX] = P[pX] - xscale;
                                                                        Qw[pY] = P[pY] - yscale;
                                                                        break;
                                                                /* 180 degrees = 06:00 */
                                                                case 4:
                                                                        Qw[pX] = P[pX];
                                                                        Qw[pY] = P[pY] - yscale;
                                                                        break;
                                                                case 3:
                                                                        Qw[pX] = P[pX] + xscale;
                                                                        Qw[pY] = P[pY] - yscale;
                                                                        break;
                                                                /* 90 degrees = 03:00 */
                                                                case 2:
                                                                        Qw[pX] = P[pX] + xscale;
                                                                        Qw[pY] = P[pY];
                                                                        break;
                                                                /* 45 degrees */
                                                                case 1:
                                                                        Qw[pX] = P[pX] + xscale;
                                                                        Qw[pY] = P[pY] + yscale;
                                                                        break;
                                                                /* 0 degrees = 00:00 */
                                                                case 0:
                                                                        Qw[pX] = P[pX];
                                                                        Qw[pY] = P[pY] + yscale;
                                                                        break;
                                                        }

                                                        /* unable to convert point to cell */
                                                        noval1 = 0;
                                                        if (rt_raster_geopoint_to_cell(
                                                                rast1,
                                                                Qw[pX], Qw[pY],
                                                                &(Qr[pX]), &(Qr[pY]),
                                                                igt1
                                                        ) != ES_NONE) {
                                                                noval1 = 1;
                                                        }
                                                        /* cell is outside bounds of grid */
                                                        else if (
                                                                (Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
                                                                (Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
                                                        ) {
                                                                noval1 = 1;
                                                        }
                                                        else if (hasnodata1 == FALSE)
                                                                val1 = 1;
                                                        /* unable to get value at cell */
                                                        else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1, &isnodata1) != ES_NONE)
                                                                noval1 = 1;

                                                        /* unable to convert point to cell */
                                                        noval2 = 0;
                                                        if (rt_raster_geopoint_to_cell(
                                                                rast2,
                                                                Qw[pX], Qw[pY],
                                                                &(Qr[pX]), &(Qr[pY]),
                                                                igt2
                                                        ) != ES_NONE) {
                                                                noval2 = 1;
                                                        }
                                                        /* cell is outside bounds of grid */
                                                        else if (
                                                                (Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
                                                                (Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
                                                        ) {
                                                                noval2 = 1;
                                                        }
                                                        else if (hasnodata2 == FALSE)
                                                                val2 = 1;
                                                        /* unable to get value at cell */
                                                        else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2, &isnodata2) != ES_NONE)
                                                                noval2 = 1;

                                                        if (!noval1) {
                                                                RASTER_DEBUGF(4, "val1 = %f", val1);
                                                        }
                                                        if (!noval2) {
                                                                RASTER_DEBUGF(4, "val2 = %f", val2);
                                                        }

                                                        /* pixels touch */
                                                        if (!noval1 && (
                                                                (hasnodata1 == FALSE) || !isnodata1
                                                        )) {
                                                                adjacent[i]++;
                                                        }
                                                        if (!noval2 && (
                                                                (hasnodata2 == FALSE) || !isnodata2
                                                        )) {
                                                                adjacent[i] += 3;
                                                        }

                                                        /* two pixel values not present */
                                                        if (noval1 || noval2) {
                                                                RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
                                                                continue;
                                                        }

                                                        /* pixels valid, so intersect */
                                                        if (
                                                                ((hasnodata1 == FALSE) || !isnodata1) &&
                                                                ((hasnodata2 == FALSE) || !isnodata2)
                                                        ) {
                                                                RASTER_DEBUG(3, "The two rasters do intersect");

                                                                return 1;
                                                        }
                                                }

                                                /* pixels touch */
                                                for (i = 0; i < 4; i++) {
                                                        RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
                                                                , i, adjacent[i], i + 4, adjacent[i + 4]);
                                                        if (adjacent[i] == 0) continue;

                                                        if (adjacent[i] + adjacent[i + 4] == 4) {
                                                                RASTER_DEBUG(3, "The two rasters touch");

                                                                return 1;
                                                        }
                                                }
                                        }
                                        else {
                                                RASTER_DEBUG(4, "outside of bounds");
                                        }
                                }
                        }
                }
        }

        return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function: