PostGIS  3.0.6dev-r@@SVN_REVISION@@

◆ 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:168
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
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:755
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
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:804
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:1221
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2235
rt_raster_get_x_scale
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
ES_NONE
@ ES_NONE
Definition: librtcore.h:180
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
rt_raster_get_y_scale
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159

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: