PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ 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.

642 {
643 int i;
644 int byHeight = 1;
645 uint32_t dimValue;
646
647 uint32_t row;
648 uint32_t rowoffset;
649 uint32_t col;
650 uint32_t coloffset;
651
652 enum line_points {X1, Y1, X2, Y2};
653 enum point {pX, pY};
654 double line1[4] = {0.};
655 double line2[4] = {0.};
656 double P[2] = {0.};
657 double Qw[2] = {0.};
658 double Qr[2] = {0.};
659 double gt1[6] = {0.};
660 double gt2[6] = {0.};
661 double igt1[6] = {0};
662 double igt2[6] = {0};
663 double d;
664 double val1;
665 int noval1;
666 int isnodata1;
667 double val2;
668 int noval2;
669 int isnodata2;
670 uint32_t adjacent[8] = {0};
671
672 double xscale;
673 double yscale;
674
675 uint16_t width1;
676 uint16_t height1;
677 uint16_t width2;
678 uint16_t height2;
679
680 width1 = rt_raster_get_width(rast1);
681 height1 = rt_raster_get_height(rast1);
682 width2 = rt_raster_get_width(rast2);
683 height2 = rt_raster_get_height(rast2);
684
685 /* sampling scale */
686 xscale = fmin(rt_raster_get_x_scale(rast1), rt_raster_get_x_scale(rast2)) / 10.;
687 yscale = fmin(rt_raster_get_y_scale(rast1), rt_raster_get_y_scale(rast2)) / 10.;
688
689 /* see if skew made rast2's rows are parallel to rast1's cols */
690 rt_raster_cell_to_geopoint(
691 rast1,
692 0, 0,
693 &(line1[X1]), &(line1[Y1]),
694 gt1
695 );
696
697 rt_raster_cell_to_geopoint(
698 rast1,
699 0, height1,
700 &(line1[X2]), &(line1[Y2]),
701 gt1
702 );
703
704 rt_raster_cell_to_geopoint(
705 rast2,
706 0, 0,
707 &(line2[X1]), &(line2[Y1]),
708 gt2
709 );
710
711 rt_raster_cell_to_geopoint(
712 rast2,
713 width2, 0,
714 &(line2[X2]), &(line2[Y2]),
715 gt2
716 );
717
718 /* Parallel lines */
719 if (FLT_EQ(((line1[X2] - line1[X1]) * (line2[Y2] - line2[Y1])),
720 ((line2[X2] - line2[X1]) * (line1[Y2] - line1[Y1]))))
721 byHeight = 0;
722
723 if (byHeight)
724 dimValue = height2;
725 else
726 dimValue = width2;
727 RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);
728
729 /* 3 x 3 search */
730 for (coloffset = 0; coloffset < 3; coloffset++) {
731 for (rowoffset = 0; rowoffset < 3; rowoffset++) {
732 /* smaller raster */
733 for (col = coloffset; col <= width1; col += 3) {
734
735 rt_raster_cell_to_geopoint(
736 rast1,
737 col, 0,
738 &(line1[X1]), &(line1[Y1]),
739 gt1
740 );
741
742 rt_raster_cell_to_geopoint(
743 rast1,
744 col, height1,
745 &(line1[X2]), &(line1[Y2]),
746 gt1
747 );
748
749 /* larger raster */
750 for (row = rowoffset; row <= dimValue; row += 3) {
751
752 if (byHeight) {
753 rt_raster_cell_to_geopoint(
754 rast2,
755 0, row,
756 &(line2[X1]), &(line2[Y1]),
757 gt2
758 );
759
760 rt_raster_cell_to_geopoint(
761 rast2,
762 width2, row,
763 &(line2[X2]), &(line2[Y2]),
764 gt2
765 );
766 }
767 else {
768 rt_raster_cell_to_geopoint(
769 rast2,
770 row, 0,
771 &(line2[X1]), &(line2[Y1]),
772 gt2
773 );
774
775 rt_raster_cell_to_geopoint(
776 rast2,
777 row, height2,
778 &(line2[X2]), &(line2[Y2]),
779 gt2
780 );
781 }
782
783 RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
784 RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
785 line1[X1], line1[Y1], line1[X2], line1[Y2]);
786 RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
787 line2[X1], line2[Y1], line2[X2], line2[Y2]);
788
789 /* intersection */
790 /* http://en.wikipedia.org/wiki/Line-line_intersection */
791 d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
792 /* no intersection */
793 if (FLT_EQ(d, 0.)) {
794 continue;
795 }
796
797 P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
798 ((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
799 P[pX] = P[pX] / d;
800
801 P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
802 ((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
803 P[pY] = P[pY] / d;
804
805 RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);
806
807 /* intersection within bounds */
808 if ((
809 (FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
810 (P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
811 ) && (
812 (FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
813 (P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
814 ) && (
815 (FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
816 (P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
817 ) && (
818 (FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
819 (P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
820 )) {
821 RASTER_DEBUG(4, "within bounds");
822
823 for (i = 0; i < 8; i++) adjacent[i] = 0;
824
825 /* test points around intersection */
826 for (i = 0; i < 8; i++) {
827 switch (i) {
828 case 7:
829 Qw[pX] = P[pX] - xscale;
830 Qw[pY] = P[pY] + yscale;
831 break;
832 /* 270 degrees = 09:00 */
833 case 6:
834 Qw[pX] = P[pX] - xscale;
835 Qw[pY] = P[pY];
836 break;
837 case 5:
838 Qw[pX] = P[pX] - xscale;
839 Qw[pY] = P[pY] - yscale;
840 break;
841 /* 180 degrees = 06:00 */
842 case 4:
843 Qw[pX] = P[pX];
844 Qw[pY] = P[pY] - yscale;
845 break;
846 case 3:
847 Qw[pX] = P[pX] + xscale;
848 Qw[pY] = P[pY] - yscale;
849 break;
850 /* 90 degrees = 03:00 */
851 case 2:
852 Qw[pX] = P[pX] + xscale;
853 Qw[pY] = P[pY];
854 break;
855 /* 45 degrees */
856 case 1:
857 Qw[pX] = P[pX] + xscale;
858 Qw[pY] = P[pY] + yscale;
859 break;
860 /* 0 degrees = 00:00 */
861 case 0:
862 Qw[pX] = P[pX];
863 Qw[pY] = P[pY] + yscale;
864 break;
865 }
866
867 /* unable to convert point to cell */
868 noval1 = 0;
869 if (rt_raster_geopoint_to_cell(
870 rast1,
871 Qw[pX], Qw[pY],
872 &(Qr[pX]), &(Qr[pY]),
873 igt1
874 ) != ES_NONE) {
875 noval1 = 1;
876 }
877 /* cell is outside bounds of grid */
878 else if ((Qr[pX] < 0 || Qr[pX] >= width1) ||
879 (Qr[pY] < 0 || Qr[pY] >= height1))
880 {
881 noval1 = 1;
882 }
883 else if (hasnodata1 == FALSE)
884 val1 = 1;
885 /* unable to get value at cell */
886 else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1, &isnodata1) != ES_NONE)
887 noval1 = 1;
888
889 /* unable to convert point to cell */
890 noval2 = 0;
891 if (rt_raster_geopoint_to_cell(
892 rast2,
893 Qw[pX], Qw[pY],
894 &(Qr[pX]), &(Qr[pY]),
895 igt2
896 ) != ES_NONE) {
897 noval2 = 1;
898 }
899 /* cell is outside bounds of grid */
900 else if ((Qr[pX] < 0 || Qr[pX] >= width2) ||
901 (Qr[pY] < 0 || Qr[pY] >= height2))
902 {
903 noval2 = 1;
904 }
905 else if (hasnodata2 == FALSE)
906 val2 = 1;
907 /* unable to get value at cell */
908 else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2, &isnodata2) != ES_NONE)
909 noval2 = 1;
910
911 if (!noval1) {
912 RASTER_DEBUGF(4, "val1 = %f", val1);
913 }
914 if (!noval2) {
915 RASTER_DEBUGF(4, "val2 = %f", val2);
916 }
917
918 /* pixels touch */
919 if (!noval1 && (
920 (hasnodata1 == FALSE) || !isnodata1
921 )) {
922 adjacent[i]++;
923 }
924 if (!noval2 && (
925 (hasnodata2 == FALSE) || !isnodata2
926 )) {
927 adjacent[i] += 3;
928 }
929
930 /* two pixel values not present */
931 if (noval1 || noval2) {
932 RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
933 continue;
934 }
935
936 /* pixels valid, so intersect */
937 if (
938 ((hasnodata1 == FALSE) || !isnodata1) &&
939 ((hasnodata2 == FALSE) || !isnodata2)
940 ) {
941 RASTER_DEBUG(3, "The two rasters do intersect");
942
943 return 1;
944 }
945 }
946
947 /* pixels touch */
948 for (i = 0; i < 4; i++) {
949 RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
950 , i, adjacent[i], i + 4, adjacent[i + 4]);
951 if (adjacent[i] == 0) continue;
952
953 if (adjacent[i] + adjacent[i + 4] == 4) {
954 RASTER_DEBUG(3, "The two rasters touch");
955
956 return 1;
957 }
958 }
959 }
960 else {
961 RASTER_DEBUG(4, "outside of bounds");
962 }
963 }
964 }
965 }
966 }
967
968 return 0;
969}
FALSE
#define FALSE
Definition dbfopen.c:72
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:304
rt_raster_cell_to_geopoint
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition rt_raster.c:637
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:308
rt_raster_geopoint_to_cell
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw, yw map point to a xr, yr raster point.
Definition rt_raster.c:686
rt_band_get_pixel
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition rt_band.c:1551
FLT_EQ
#define FLT_EQ(x, y)
Definition librtcore.h:2436
rt_raster_get_x_scale
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition rt_raster.c:154
ES_NONE
@ ES_NONE
Definition librtcore.h:182
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition rt_raster.c:133
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition rt_raster.c:125
rt_raster_get_y_scale
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition rt_raster.c:163

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().

Here is the call graph for this function:
Here is the caller graph for this function: