verkkoalgoritmi

Ari Jolma ajolma at water.hut.fi
Fri Nov 19 03:50:25 CST 1999


Tämä ei ole varsinaisesti Perl asia vaan yleisempi mutta silti. Puuhailen
Perl-moduulin
kimpussa, jonka on tarkoitus olla a) yleinen rastereiden (tai mitä termiä
haluaakin käyttää,
kyse on matriiseista) käpistelytyökalu, b) yleinen GIS työkalu ja c)
yleinen hydrologinen
GIS työkalu. Nämä rasterit ovat isoja (esim. 3000x3000) ja siksi algoritmit
on pääasiassa C:llä ja Perlillä on tehty näppärä oliopohjainen ja
operaattorit määrittelevä API.

Mutta siis siihen algoritmiin, joka tuottaa ongelmia koska ne
implementaatiot jotka olen
saanut aikaan ovat hitaita - luokkaa >10 tuntia isolle rasterille ajettuna
eli esim. 10 x nopeutus auttaisi tosi paljon (ajelen näitä omilla
PIII:lla). Mulla on siis FDG (Flow Direction Grid), jonka jokaisessa
pisteessä (i,j, jotka kuten matriiseissa) on kokonaisluku 0..9, 0
tarkoittaa, että solu on "pit" eli kuoppa, 1 tarkoittaa, että vesi virtaa
tästä pohjoiseen eli soluun (i-1,j), 2 tarkoittaa, että vesi virtaa
koilliseen eli soluun (i-1,j+1), jne. Ideana olisi poistaa painanteet eli
alueet, joissa vesi valuu soluun, joka ei ole gridin reunalla. Olen tehnyt 
sen, niin, että (i) etsin alueen reunalta alimman solun (ii) muutan polun
siitä pitiin
kulkemaan toisinpäin ja (iii) sen itsensä purkautumaan alueen ulkopuolelle
alimpaan
naapuriinsa.

Olen sillä tavalla hölmö, että minusta näitä on hauska väsätä itse ja tällä
kertaa se voi olla
syynä siihen että tulos on noin hidas. Ehkäpä algoritmikirjallisuudessa
olisi parempi ratkaisu, ja joku voisi viitata minut siihen tai antaa
mielestään nopean koodin niin testaan.
Nopea ratkaisu voisi olla sellainen, joka kävisi koko alueen läpi ja hakisi
reunapisteet
siinä samalla.

Vielä sen verran että FDG lasketaan DEMistä, joka on korkeusmalli mutta
tällä seikalla ei ole tämän algoritmin kanssa tekemistä.

Toivottavasti en pilaa nyt kenenkään viikonloppua ;-)

Ari

tässä mun tämänhetkinen versio:

int fixpit(GRID *fdg, GRID *dem, GRID *mark, int i, int j) {
  int i0,j0,i2,j2,ilbp,jlbp,ilbpn,jlbpn,d;
  REAL lbp,lbpn;
#ifdef DEBUG_ENABLED
  if (debug) fprintf(stderr,"fix pit at %i,%i\n",i,j);
#endif
  memset(mark->sdata,0,mark->M*mark->N*sizeof(INTEGER));
  markbasin(fdg, mark, i, j);
#ifdef DEBUG_ENABLED
  if (debug>1) fprintf(stderr,"basin marked %i,%i\n",i,j);

#endif
  /* walk the border, start from top and walk clockwise */
  i2 = i;
  j2 = j;
  while (i2 > 0 AND sgd(mark,i2-1,j2)) {
#ifdef DEBUG_ENABLED
    if (debug>3) fprintf(stderr,"up to border, at cell %i,%i\n",i2,j2);
#endif
    i2--;
  }
  /* d is the direction where we go */
  d = 3;
#ifdef DEBUG_ENABLED
  if (debug>3) fprintf(stderr,"start walking, at cell %i,%i, going
%i\n",i2,j2,d);
#endif
  ilbp = i0 = i2;
  jlbp = j0 = j2;
  lbp = fgd(dem,i2,j2);
  do {
    /* test clockwise, the first dir to test depends on our current
direction */
    /* test direction is d2 */
    int d2 = -1;
    int done = -1;
#ifdef DEBUG_ENABLED
    if (debug>3) fprintf(stderr,"walking, at cell %i,%i, going %i, current
low is %i,%i %f\n",i2,j2,d,ilbp,jlbp,lbp);
#endif
    switch (d) {
    sagain:
    case 1:
    case 2:
      {
	if (i2 > 0 AND j2 > 0 AND sgd(mark,i2-1,j2-1)) {
	  d2 = 8;
	  i2--;
	  j2--;
	  break;
	}
	if (i2 > 0 AND sgd(mark,i2-1,j2)) {
	  d2 = 1;
	  i2--;
	  break;
	}
      }
    case 3:
    case 4:
      {
	if (i2 > 0 AND j2 < mark->N-1 AND sgd(mark,i2-1,j2+1)) {
	  d2 = 2;
	  i2--;
	  j2++;
	  break;
	}
	if (j2 < mark->N-1 AND sgd(mark,i2,j2+1)) {
	  d2 = 3;
	  j2++;
	  break;
	}
      }
    case 5:
    case 6:
      {
	if (i2 < mark->M-1 AND j2 < mark->N-1 AND sgd(mark,i2+1,j2+1)) {
	  d2 = 4;
	  i2++;
	  j2++;
	  break;
	}
	if (i2 < mark->M-1 AND sgd(mark,i2+1,j2)) {
	  d2 = 5;
	  i2++;
	  break;
	}
      }
    case 7:
    case 8:
      {
	if (i2 < mark->M-1 AND j2 > 0 AND sgd(mark,i2+1,j2-1)) {
	  d2 = 6;
	  i2++;
	  j2--;
	  break;
	}
	if (j2 > 0 AND sgd(mark,i2,j2-1)) {
	  d2 = 7;
	  j2--;
	  break;
	}
      }
    }
    /* if new direction is not found, try once more */
    done++;
#ifdef DEBUG_ENABLED
    if (debug>3) fprintf(stderr,"do I have the new dir (%i)? %i\n",d2,done);
#endif
    if (d2 == -1 AND !done) goto sagain;
    d = d2;
    if (fgd(dem,i2,j2) < lbp) {
      lbp = fgd(dem,i2,j2);
      ilbp = i2;
      jlbp = j2;
    }

  } while (!(i2 == i0 AND j2 == j0));
#ifdef DEBUG_ENABLED
  if (debug>1) fprintf(stderr,"lowest border point is at %i,%i\n",ilbp,jlbp);
#endif
  /* 
     if the lowest border point is at the border of the whole grid then 
     that will become a pit, else search for the lowest neighbor not in the
     depression
  */
  ilbpn = -1;
  jlbpn = -1;
  if (ilbp > 0 AND ilbp < mark->M-1 AND jlbp > 0 AND jlbp < mark->N-1) {
    lbpn = 99999999.9;
    i2 = ilbp;
    j2 = jlbp;
    i2--;
    if (i2 >= 0 AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) < lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    j2++;
    if (i2 >= 0 AND j2 < fdg->N AND !sgd(mark,i2,j2) AND fgd(dem,i2,j2) <
lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    i2++;
    if (j2 < fdg->N AND !sgd(mark,i2,j2)  AND fgd(dem,i2,j2) < lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    i2++;
    if (i2 < fdg->M AND j2 < fdg->N AND !sgd(mark,i2,j2)  AND
fgd(dem,i2,j2) < lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    j2--;
    if (i2 < fdg->M AND !sgd(mark,i2,j2)  AND fgd(dem,i2,j2) < lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    j2--;
    if (i2 < fdg->M AND j2 >= 0 AND !sgd(mark,i2,j2)  AND fgd(dem,i2,j2) <
lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    i2--;
    if (j2 >= 0 AND !sgd(mark,i2,j2)  AND fgd(dem,i2,j2) < lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    i2--;
    if (i2 >= 0 AND j2 >= 0 AND !sgd(mark,i2,j2)  AND fgd(dem,i2,j2) < lbpn) {
      lbpn = fgd(dem,i2,j2);
      ilbpn = i2;
      jlbpn = j2;
    }
    if (ilbpn == -1) {
      fprintf(stderr,"fixpit: FAILED: the lowest border point of a
depression at point (%i,%i) does not have a neigbor outside the
depression\n",i,j);
      return 0;
    }
#ifdef DEBUG_ENABLED
    if (debug>1) fprintf(stderr,"escape from the depression:
%i,%i,%i,%i\n",ilbp,jlbp,ilbpn,jlbpn);
#endif
  }
  if (ilbp != i OR jlbp != j) {
    /* invert the path from the lowest border point to the pit point */
    int im = ilbp;
    int jm = jlbp;
    int da = 0;
    int db = sgd(fdg,im,jm);
    int tmp;
    i2 = ilbp;
    j2 = jlbp;
    do {
      movetodir(i2, j2, db);
      inversedir(da, tmp);
      sgd(fdg,im,jm) = tmp;
      im = i2;
      jm = j2;
      da = db;
      db = sgd(fdg,im,jm);
    } while (db);
    inversedir(da, tmp);
    sgd(fdg,im,jm) = tmp;
  }
  
  /* a new pit point at the border of the grid or escape from the
depression */
  if (ilbpn == -1) sgd(fdg,ilbp,jlbp) = 0;
  else {
    if (jlbp == jlbpn) {
      if (ilbpn-ilbp == 1) sgd(fdg,ilbp,jlbp) = 5;
      else if (ilbp-ilbpn == 1) sgd(fdg,ilbp,jlbp) = 1;
    } else if (jlbp-jlbpn == 1) {
      if (ilbpn == ilbp) sgd(fdg,ilbp,jlbp) = 7;
      else if (ilbpn-ilbp == 1) sgd(fdg,ilbp,jlbp) = 6;
      else if (ilbp-ilbpn == 1) sgd(fdg,ilbp,jlbp) = 8;
    } else if (jlbpn-jlbp == 1) {
      if (ilbpn == ilbp) sgd(fdg,ilbp,jlbp) = 3;
      else if (ilbpn-ilbp == 1) sgd(fdg,ilbp,jlbp) = 4;
      else if (ilbp-ilbpn == 1) sgd(fdg,ilbp,jlbp) = 2;
    }
  }
#ifdef DEBUG_ENABLED
    if (debug>1) fprintf(stderr,"pit fixed, hopefully\n");
#endif
  return 1;
}





More information about the Helsinki-pm mailing list