PostGIS  2.5.1dev-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.

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

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 */
691  rast1,
692  0, 0,
693  &(line1[X1]), &(line1[Y1]),
694  gt1
695  );
696 
698  rast1,
699  0, height1,
700  &(line1[X2]), &(line1[Y2]),
701  gt1
702  );
703 
705  rast2,
706  0, 0,
707  &(line2[X1]), &(line2[Y1]),
708  gt2
709  );
710 
712  rast2,
713  width2, 0,
714  &(line2[X2]), &(line2[Y2]),
715  gt2
716  );
717 
718  /* parallel vertically */
719  if (FLT_EQ(line1[X2] - line1[X1], 0.) && FLT_EQ(line2[X2] - line2[X1], 0.))
720  byHeight = 0;
721  /* parallel */
722  else if (FLT_EQ(((line1[Y2] - line1[Y1]) / (line1[X2] - line1[X1])), ((line2[Y2] - line2[Y1]) / (line2[X2] - line2[X1]))))
723  byHeight = 0;
724 
725  if (byHeight)
726  dimValue = height2;
727  else
728  dimValue = width2;
729  RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);
730 
731  /* 3 x 3 search */
732  for (coloffset = 0; coloffset < 3; coloffset++) {
733  for (rowoffset = 0; rowoffset < 3; rowoffset++) {
734  /* smaller raster */
735  for (col = coloffset; col <= width1; col += 3) {
736 
738  rast1,
739  col, 0,
740  &(line1[X1]), &(line1[Y1]),
741  gt1
742  );
743 
745  rast1,
746  col, height1,
747  &(line1[X2]), &(line1[Y2]),
748  gt1
749  );
750 
751  /* larger raster */
752  for (row = rowoffset; row <= dimValue; row += 3) {
753 
754  if (byHeight) {
756  rast2,
757  0, row,
758  &(line2[X1]), &(line2[Y1]),
759  gt2
760  );
761 
763  rast2,
764  width2, row,
765  &(line2[X2]), &(line2[Y2]),
766  gt2
767  );
768  }
769  else {
771  rast2,
772  row, 0,
773  &(line2[X1]), &(line2[Y1]),
774  gt2
775  );
776 
778  rast2,
779  row, height2,
780  &(line2[X2]), &(line2[Y2]),
781  gt2
782  );
783  }
784 
785  RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
786  RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
787  line1[X1], line1[Y1], line1[X2], line1[Y2]);
788  RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
789  line2[X1], line2[Y1], line2[X2], line2[Y2]);
790 
791  /* intersection */
792  /* http://en.wikipedia.org/wiki/Line-line_intersection */
793  d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
794  /* no intersection */
795  if (FLT_EQ(d, 0.)) {
796  continue;
797  }
798 
799  P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
800  ((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
801  P[pX] = P[pX] / d;
802 
803  P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
804  ((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
805  P[pY] = P[pY] / d;
806 
807  RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);
808 
809  /* intersection within bounds */
810  if ((
811  (FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
812  (P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
813  ) && (
814  (FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
815  (P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
816  ) && (
817  (FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
818  (P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
819  ) && (
820  (FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
821  (P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
822  )) {
823  RASTER_DEBUG(4, "within bounds");
824 
825  for (i = 0; i < 8; i++) adjacent[i] = 0;
826 
827  /* test points around intersection */
828  for (i = 0; i < 8; i++) {
829  switch (i) {
830  case 7:
831  Qw[pX] = P[pX] - xscale;
832  Qw[pY] = P[pY] + yscale;
833  break;
834  /* 270 degrees = 09:00 */
835  case 6:
836  Qw[pX] = P[pX] - xscale;
837  Qw[pY] = P[pY];
838  break;
839  case 5:
840  Qw[pX] = P[pX] - xscale;
841  Qw[pY] = P[pY] - yscale;
842  break;
843  /* 180 degrees = 06:00 */
844  case 4:
845  Qw[pX] = P[pX];
846  Qw[pY] = P[pY] - yscale;
847  break;
848  case 3:
849  Qw[pX] = P[pX] + xscale;
850  Qw[pY] = P[pY] - yscale;
851  break;
852  /* 90 degrees = 03:00 */
853  case 2:
854  Qw[pX] = P[pX] + xscale;
855  Qw[pY] = P[pY];
856  break;
857  /* 45 degrees */
858  case 1:
859  Qw[pX] = P[pX] + xscale;
860  Qw[pY] = P[pY] + yscale;
861  break;
862  /* 0 degrees = 00:00 */
863  case 0:
864  Qw[pX] = P[pX];
865  Qw[pY] = P[pY] + yscale;
866  break;
867  }
868 
869  /* unable to convert point to cell */
870  noval1 = 0;
872  rast1,
873  Qw[pX], Qw[pY],
874  &(Qr[pX]), &(Qr[pY]),
875  igt1
876  ) != ES_NONE) {
877  noval1 = 1;
878  }
879  /* cell is outside bounds of grid */
880  else if (
881  (Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
882  (Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
883  ) {
884  noval1 = 1;
885  }
886  else if (hasnodata1 == FALSE)
887  val1 = 1;
888  /* unable to get value at cell */
889  else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1, &isnodata1) != ES_NONE)
890  noval1 = 1;
891 
892  /* unable to convert point to cell */
893  noval2 = 0;
895  rast2,
896  Qw[pX], Qw[pY],
897  &(Qr[pX]), &(Qr[pY]),
898  igt2
899  ) != ES_NONE) {
900  noval2 = 1;
901  }
902  /* cell is outside bounds of grid */
903  else if (
904  (Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
905  (Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
906  ) {
907  noval2 = 1;
908  }
909  else if (hasnodata2 == FALSE)
910  val2 = 1;
911  /* unable to get value at cell */
912  else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2, &isnodata2) != ES_NONE)
913  noval2 = 1;
914 
915  if (!noval1) {
916  RASTER_DEBUGF(4, "val1 = %f", val1);
917  }
918  if (!noval2) {
919  RASTER_DEBUGF(4, "val2 = %f", val2);
920  }
921 
922  /* pixels touch */
923  if (!noval1 && (
924  (hasnodata1 == FALSE) || !isnodata1
925  )) {
926  adjacent[i]++;
927  }
928  if (!noval2 && (
929  (hasnodata2 == FALSE) || !isnodata2
930  )) {
931  adjacent[i] += 3;
932  }
933 
934  /* two pixel values not present */
935  if (noval1 || noval2) {
936  RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
937  continue;
938  }
939 
940  /* pixels valid, so intersect */
941  if (
942  ((hasnodata1 == FALSE) || !isnodata1) &&
943  ((hasnodata2 == FALSE) || !isnodata2)
944  ) {
945  RASTER_DEBUG(3, "The two rasters do intersect");
946 
947  return 1;
948  }
949  }
950 
951  /* pixels touch */
952  for (i = 0; i < 4; i++) {
953  RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
954  , i, adjacent[i], i + 4, adjacent[i + 4]);
955  if (adjacent[i] == 0) continue;
956 
957  if (adjacent[i] + adjacent[i + 4] == 4) {
958  RASTER_DEBUG(3, "The two rasters touch");
959 
960  return 1;
961  }
962  }
963  }
964  else {
965  RASTER_DEBUG(4, "outside of bounds");
966  }
967  }
968  }
969  }
970  }
971 
972  return 0;
973 }
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
unsigned int uint32_t
Definition: uthash.h:78
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
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
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
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:806
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
Here is the call graph for this function:
Here is the caller graph for this function: