PostGIS  2.5.7dev-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 */
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 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 
736  rast1,
737  col, 0,
738  &(line1[X1]), &(line1[Y1]),
739  gt1
740  );
741 
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) {
754  rast2,
755  0, row,
756  &(line2[X1]), &(line2[Y1]),
757  gt2
758  );
759 
761  rast2,
762  width2, row,
763  &(line2[X2]), &(line2[Y2]),
764  gt2
765  );
766  }
767  else {
769  rast2,
770  row, 0,
771  &(line2[X1]), &(line2[Y1]),
772  gt2
773  );
774 
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;
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 (
879  (Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
880  (Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
881  ) {
882  noval1 = 1;
883  }
884  else if (hasnodata1 == FALSE)
885  val1 = 1;
886  /* unable to get value at cell */
887  else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1, &isnodata1) != ES_NONE)
888  noval1 = 1;
889 
890  /* unable to convert point to cell */
891  noval2 = 0;
893  rast2,
894  Qw[pX], Qw[pY],
895  &(Qr[pX]), &(Qr[pY]),
896  igt2
897  ) != ES_NONE) {
898  noval2 = 1;
899  }
900  /* cell is outside bounds of grid */
901  else if (
902  (Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
903  (Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
904  ) {
905  noval2 = 1;
906  }
907  else if (hasnodata2 == FALSE)
908  val2 = 1;
909  /* unable to get value at cell */
910  else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2, &isnodata2) != ES_NONE)
911  noval2 = 1;
912 
913  if (!noval1) {
914  RASTER_DEBUGF(4, "val1 = %f", val1);
915  }
916  if (!noval2) {
917  RASTER_DEBUGF(4, "val2 = %f", val2);
918  }
919 
920  /* pixels touch */
921  if (!noval1 && (
922  (hasnodata1 == FALSE) || !isnodata1
923  )) {
924  adjacent[i]++;
925  }
926  if (!noval2 && (
927  (hasnodata2 == FALSE) || !isnodata2
928  )) {
929  adjacent[i] += 3;
930  }
931 
932  /* two pixel values not present */
933  if (noval1 || noval2) {
934  RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
935  continue;
936  }
937 
938  /* pixels valid, so intersect */
939  if (
940  ((hasnodata1 == FALSE) || !isnodata1) &&
941  ((hasnodata2 == FALSE) || !isnodata2)
942  ) {
943  RASTER_DEBUG(3, "The two rasters do intersect");
944 
945  return 1;
946  }
947  }
948 
949  /* pixels touch */
950  for (i = 0; i < 4; i++) {
951  RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
952  , i, adjacent[i], i + 4, adjacent[i + 4]);
953  if (adjacent[i] == 0) continue;
954 
955  if (adjacent[i] + adjacent[i + 4] == 4) {
956  RASTER_DEBUG(3, "The two rasters touch");
957 
958  return 1;
959  }
960  }
961  }
962  else {
963  RASTER_DEBUG(4, "outside of bounds");
964  }
965  }
966  }
967  }
968  }
969 
970  return 0;
971 }
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
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
#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
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
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
@ ES_NONE
Definition: librtcore.h:180
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
unsigned int uint32_t
Definition: uthash.h:78

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: