PostGIS  2.1.10dev-r@@SVN_REVISION@@
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 11598 of file rt_api.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().

11603  {
11604  int i;
11605  int byHeight = 1;
11606  uint32_t dimValue;
11607 
11608  uint32_t row;
11609  uint32_t rowoffset;
11610  uint32_t col;
11611  uint32_t coloffset;
11612 
11613  enum line_points {X1, Y1, X2, Y2};
11614  enum point {pX, pY};
11615  double line1[4] = {0.};
11616  double line2[4] = {0.};
11617  double P[2] = {0.};
11618  double Qw[2] = {0.};
11619  double Qr[2] = {0.};
11620  double gt1[6] = {0.};
11621  double gt2[6] = {0.};
11622  double igt1[6] = {0};
11623  double igt2[6] = {0};
11624  double d;
11625  double val1;
11626  int noval1;
11627  int isnodata1;
11628  double val2;
11629  int noval2;
11630  int isnodata2;
11631  uint32_t adjacent[8] = {0};
11632 
11633  double xscale;
11634  double yscale;
11635 
11636  uint16_t width1;
11637  uint16_t height1;
11638  uint16_t width2;
11639  uint16_t height2;
11640 
11641  width1 = rt_raster_get_width(rast1);
11642  height1 = rt_raster_get_height(rast1);
11643  width2 = rt_raster_get_width(rast2);
11644  height2 = rt_raster_get_height(rast2);
11645 
11646  /* sampling scale */
11647  xscale = fmin(rt_raster_get_x_scale(rast1), rt_raster_get_x_scale(rast2)) / 10.;
11648  yscale = fmin(rt_raster_get_y_scale(rast1), rt_raster_get_y_scale(rast2)) / 10.;
11649 
11650  /* see if skew made rast2's rows are parallel to rast1's cols */
11652  rast1,
11653  0, 0,
11654  &(line1[X1]), &(line1[Y1]),
11655  gt1
11656  );
11657 
11659  rast1,
11660  0, height1,
11661  &(line1[X2]), &(line1[Y2]),
11662  gt1
11663  );
11664 
11666  rast2,
11667  0, 0,
11668  &(line2[X1]), &(line2[Y1]),
11669  gt2
11670  );
11671 
11673  rast2,
11674  width2, 0,
11675  &(line2[X2]), &(line2[Y2]),
11676  gt2
11677  );
11678 
11679  /* parallel vertically */
11680  if (FLT_EQ(line1[X2] - line1[X1], 0.) && FLT_EQ(line2[X2] - line2[X1], 0.))
11681  byHeight = 0;
11682  /* parallel */
11683  else if (FLT_EQ(((line1[Y2] - line1[Y1]) / (line1[X2] - line1[X1])), ((line2[Y2] - line2[Y1]) / (line2[X2] - line2[X1]))))
11684  byHeight = 0;
11685 
11686  if (byHeight)
11687  dimValue = height2;
11688  else
11689  dimValue = width2;
11690  RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);
11691 
11692  /* 3 x 3 search */
11693  for (coloffset = 0; coloffset < 3; coloffset++) {
11694  for (rowoffset = 0; rowoffset < 3; rowoffset++) {
11695  /* smaller raster */
11696  for (col = coloffset; col <= width1; col += 3) {
11697 
11699  rast1,
11700  col, 0,
11701  &(line1[X1]), &(line1[Y1]),
11702  gt1
11703  );
11704 
11706  rast1,
11707  col, height1,
11708  &(line1[X2]), &(line1[Y2]),
11709  gt1
11710  );
11711 
11712  /* larger raster */
11713  for (row = rowoffset; row <= dimValue; row += 3) {
11714 
11715  if (byHeight) {
11717  rast2,
11718  0, row,
11719  &(line2[X1]), &(line2[Y1]),
11720  gt2
11721  );
11722 
11724  rast2,
11725  width2, row,
11726  &(line2[X2]), &(line2[Y2]),
11727  gt2
11728  );
11729  }
11730  else {
11732  rast2,
11733  row, 0,
11734  &(line2[X1]), &(line2[Y1]),
11735  gt2
11736  );
11737 
11739  rast2,
11740  row, height2,
11741  &(line2[X2]), &(line2[Y2]),
11742  gt2
11743  );
11744  }
11745 
11746  RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
11747  RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
11748  line1[X1], line1[Y1], line1[X2], line1[Y2]);
11749  RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
11750  line2[X1], line2[Y1], line2[X2], line2[Y2]);
11751 
11752  /* intersection */
11753  /* http://en.wikipedia.org/wiki/Line-line_intersection */
11754  d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
11755  /* no intersection */
11756  if (FLT_EQ(d, 0.)) {
11757  continue;
11758  }
11759 
11760  P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
11761  ((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
11762  P[pX] = P[pX] / d;
11763 
11764  P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
11765  ((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
11766  P[pY] = P[pY] / d;
11767 
11768  RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);
11769 
11770  /* intersection within bounds */
11771  if ((
11772  (FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
11773  (P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
11774  ) && (
11775  (FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
11776  (P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
11777  ) && (
11778  (FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
11779  (P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
11780  ) && (
11781  (FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
11782  (P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
11783  )) {
11784  RASTER_DEBUG(4, "within bounds");
11785 
11786  for (i = 0; i < 8; i++) adjacent[i] = 0;
11787 
11788  /* test points around intersection */
11789  for (i = 0; i < 8; i++) {
11790  switch (i) {
11791  case 7:
11792  Qw[pX] = P[pX] - xscale;
11793  Qw[pY] = P[pY] + yscale;
11794  break;
11795  /* 270 degrees = 09:00 */
11796  case 6:
11797  Qw[pX] = P[pX] - xscale;
11798  Qw[pY] = P[pY];
11799  break;
11800  case 5:
11801  Qw[pX] = P[pX] - xscale;
11802  Qw[pY] = P[pY] - yscale;
11803  break;
11804  /* 180 degrees = 06:00 */
11805  case 4:
11806  Qw[pX] = P[pX];
11807  Qw[pY] = P[pY] - yscale;
11808  break;
11809  case 3:
11810  Qw[pX] = P[pX] + xscale;
11811  Qw[pY] = P[pY] - yscale;
11812  break;
11813  /* 90 degrees = 03:00 */
11814  case 2:
11815  Qw[pX] = P[pX] + xscale;
11816  Qw[pY] = P[pY];
11817  break;
11818  /* 45 degrees */
11819  case 1:
11820  Qw[pX] = P[pX] + xscale;
11821  Qw[pY] = P[pY] + yscale;
11822  break;
11823  /* 0 degrees = 00:00 */
11824  case 0:
11825  Qw[pX] = P[pX];
11826  Qw[pY] = P[pY] + yscale;
11827  break;
11828  }
11829 
11830  /* unable to convert point to cell */
11831  noval1 = 0;
11833  rast1,
11834  Qw[pX], Qw[pY],
11835  &(Qr[pX]), &(Qr[pY]),
11836  igt1
11837  ) != ES_NONE) {
11838  noval1 = 1;
11839  }
11840  /* cell is outside bounds of grid */
11841  else if (
11842  (Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
11843  (Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
11844  ) {
11845  noval1 = 1;
11846  }
11847  else if (hasnodata1 == FALSE)
11848  val1 = 1;
11849  /* unable to get value at cell */
11850  else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1, &isnodata1) != ES_NONE)
11851  noval1 = 1;
11852 
11853  /* unable to convert point to cell */
11854  noval2 = 0;
11856  rast2,
11857  Qw[pX], Qw[pY],
11858  &(Qr[pX]), &(Qr[pY]),
11859  igt2
11860  ) != ES_NONE) {
11861  noval2 = 1;
11862  }
11863  /* cell is outside bounds of grid */
11864  else if (
11865  (Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
11866  (Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
11867  ) {
11868  noval2 = 1;
11869  }
11870  else if (hasnodata2 == FALSE)
11871  val2 = 1;
11872  /* unable to get value at cell */
11873  else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2, &isnodata2) != ES_NONE)
11874  noval2 = 1;
11875 
11876  if (!noval1) {
11877  RASTER_DEBUGF(4, "val1 = %f", val1);
11878  }
11879  if (!noval2) {
11880  RASTER_DEBUGF(4, "val2 = %f", val2);
11881  }
11882 
11883  /* pixels touch */
11884  if (!noval1 && (
11885  (hasnodata1 == FALSE) || !isnodata1
11886  )) {
11887  adjacent[i]++;
11888  }
11889  if (!noval2 && (
11890  (hasnodata2 == FALSE) || !isnodata2
11891  )) {
11892  adjacent[i] += 3;
11893  }
11894 
11895  /* two pixel values not present */
11896  if (noval1 || noval2) {
11897  RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
11898  continue;
11899  }
11900 
11901  /* pixels valid, so intersect */
11902  if (
11903  ((hasnodata1 == FALSE) || !isnodata1) &&
11904  ((hasnodata2 == FALSE) || !isnodata2)
11905  ) {
11906  RASTER_DEBUG(3, "The two rasters do intersect");
11907 
11908  return 1;
11909  }
11910  }
11911 
11912  /* pixels touch */
11913  for (i = 0; i < 4; i++) {
11914  RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
11915  , i, adjacent[i], i + 4, adjacent[i + 4]);
11916  if (adjacent[i] == 0) continue;
11917 
11918  if (adjacent[i] + adjacent[i + 4] == 4) {
11919  RASTER_DEBUG(3, "The two rasters touch");
11920 
11921  return 1;
11922  }
11923  }
11924  }
11925  else {
11926  RASTER_DEBUG(4, "outside of bounds");
11927  }
11928  }
11929  }
11930  }
11931  }
11932 
11933  return 0;
11934 }
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_api.c:6054
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_api.c:5455
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_api.c:5434
#define FALSE
Definition: dbfopen.c:169
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_api.c:5464
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_api.c:6105
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_api.c:5426

Here is the call graph for this function:

Here is the caller graph for this function: