Blockstructured Adaptive Mesh Refinement in object-oriented C++

---

       

---

Cluster3.C

Go to the documentation of this file.

#include "DAGH.h"
#include "values.h"

#ifdef __KCC
 template class GridData(2)<short>;
 template class GridData(3)<short>;
 template class GridData(1)<int>;
#endif

#define MAXPOINTS                    (10000)
#define MAXPOINTS2               (100000000)
#define MIN_EFFICIENCY                 (0.4)
#define MINOFF                           (4)
#define ITHRES                           (2)

void Cluster_Prune(GridData(1)<int> &sig, int& pl, int& pu, int bw) {
  BeginFastIndex1(sig, sig.bbox(), sig.data(), int);
  pl = sig.bbox().lower(0);
  pu = sig.bbox().upper(0);
  if (pu - pl < bw) {
    return;
  }
  while (FastIndex1(sig,pl) == 0 && pl < sig.bbox().upper(0)) 
    pl += sig.bbox().stepsize(0);
  while (FastIndex1(sig,pu) == 0 && pu > sig.bbox().lower(0)) 
    pu -= sig.bbox().stepsize(0);
  if (pu+sig.bbox().stepsize(0)-pl < bw) {
    int diff = bw - (pu+sig.bbox().stepsize(0)-pl);
    int addl, addu;
    int udist = sig.bbox().upper(0) - pu;
    int ldist = pl - sig.bbox().lower(0);

    addl = ((diff/sig.bbox().stepsize(0))/2) * sig.bbox().stepsize(0);
    addu = diff - addl;
    if (addu > udist) {
      addu = udist;
      addl = diff - addu;
    }
    if (addl > ldist) {
      addl = ldist;
      addu = diff - addu;
    }

    pu += addu;
    pl -= addl;
  }
  EndFastIndex1(sig);
}

void Cluster_Slice(GridData(1)<int> &sig, int& slice, int& str, int bw) {
  int i, l, u, s, t, a;         // ctr, low, up, step, str, answer
  l = sig.bbox().lower(0);
  u = sig.bbox().upper(0);
  s = sig.bbox().stepsize(0);

  /* Since this box is already pruned, we don't want to split it if
     it is less than 2 the min width */
  if (u+s-l < 2*bw) {
    str   = -1;
    slice = -1;
    return;
  }

  GridData(1)<int> lap(l,u,s);

  BeginFastIndex1(sig, sig.bbox(), sig.data(), int);
  BeginFastIndex1(lap, lap.bbox(), lap.data(), int);

  t = -1;
  a = -1;
  
  /* Look for a zero near the center */
  for (i=l+bw-s;i<=u-bw;i+=s) {
    if (FastIndex1(sig,i) == 0) {
      int di = (i-l < u-i ? i-l : u-i);
      if (di > t) {
        a = i;
        t = di;
      }
    }
  }
  if (a != -1) {
    slice = a;
    str = 2*MAXPOINTS2+t;
    return;
  }

  /* Evaluate the laplacian etc */
  int malap = 0;
  for (i=l+s; i<=u-s; i+=s) {
    lap(i) = FastIndex1(sig,i+s) + FastIndex1(sig,i-s) - 2*FastIndex1(sig,i);
    int alap = (FastIndex1(lap,i)>0 ? FastIndex1(lap,i) : -FastIndex1(lap,i));
    if (alap > malap) malap = alap;
  }

  t = -1;
  a = -1;

  /* And find the strongest second derivative */
  int tts = -1;
  if (malap > 0) {
    int ms = Max(MINOFF*s,bw);
    for (i=l+ms-s; i<=u-ms; i+=s)
      if (FastIndex1(lap,i-s)*FastIndex1(lap,i) < 0) {  
        int ts = ((FastIndex1(lap,i-s)-FastIndex1(lap,i))>0 ? 
                  FastIndex1(lap,i-s)-FastIndex1(lap,i) : 
                  -(FastIndex1(lap,i-s)-FastIndex1(lap,i)));
        int di = (i-l < u-i ? i-l : u-i);
        if (ts==tts && di>t) {
          a = i;
          t = di;
        }
        if (ts>tts) {
          tts = ts;
          a = i;
          t = di;
        }
      }
    if (a!=-1 && tts>=ITHRES) 
      t = MAXPOINTS2+tts*MAXPOINTS+t;
  }

  /* Use bisection if no slice has been found */
  if (a==-1 || tts<ITHRES) {
    t = (sig.size()/2)*s-s;
    a = l+t;
  }
  EndFastIndex1(sig);
  EndFastIndex1(lap);

  slice = a;
  str = t;
}

void Cluster3R(GridData(3)<short> &flag, 
               const BBox& flagbbox,
               double Efficiency,
               int MinWidth,
               int MaxWidth,
               BBoxList &recurseOnThis) {
  BBoxList result, tmpresult;
  int l[3], u[3], s[3];
  int i;
  for (i=0;i<3;i++) s[i] = flagbbox.stepsize(i);
  if (MaxWidth<2*MinWidth) MaxWidth=MAXINT;

  //cout << "Clustering " << flagbbox <<"\n";
  //cout.flush();

  /* Declare and fill in Signature Planes in each direction */
  Array(1)<int> i_sig(flagbbox.lower(0), 
                      flagbbox.upper(0),
                      flagbbox.stepsize(0));
  Array(1)<int> j_sig(flagbbox.lower(1), 
                      flagbbox.upper(1), 
                      flagbbox.stepsize(1));
  Array(1)<int> k_sig(flagbbox.lower(2), 
                      flagbbox.upper(2),
                      flagbbox.stepsize(2));
  i_sig = 0;
  j_sig = 0;
  k_sig = 0;
  BeginFastIndex3(flag, flag.bbox(), flag.data(), short);
  BeginFastIndex1(i_sig, i_sig.bbox(), i_sig.data(), int);
  BeginFastIndex1(j_sig, j_sig.bbox(), j_sig.data(), int);
  BeginFastIndex1(k_sig, k_sig.bbox(), k_sig.data(), int);

  int flagged = 0;
  for_3(ii, jj, kk, flagbbox, flagbbox.stepsize())
    if (FastIndex3(flag,ii,jj,kk)>0) {
      FastIndex1(i_sig,ii)++;
      FastIndex1(j_sig,jj)++;
      FastIndex1(k_sig,kk)++;
      flagged++;
    }
  end_for;

  /* Check for pruning */
  int ipl, ipu, jpl, jpu, kpl, kpu;
  Cluster_Prune(i_sig,ipl,ipu,MinWidth*s[0]);
  Cluster_Prune(j_sig,jpl,jpu,MinWidth*s[1]);
  Cluster_Prune(k_sig,kpl,kpu,MinWidth*s[2]);

  Coords pl(3,ipl,jpl,kpl);
  Coords pu(3,ipu,jpu,kpu);

  if (! ((pl == flagbbox.lower()) &&
         (pu == flagbbox.upper()))) {
    Coords ps = Coords(3,s);
    BBox newB(pl,pu,ps);
    // cout <<"  PRUNE :" << newB <<"\n";
    // cout.flush();
    /* Check my efficiency! */
    flagged = 0;
    for_3(ii, jj, kk, newB, newB.stepsize())
      if (FastIndex3(flag,ii,jj,kk)>0) flagged++;
    end_for
    if (((double)flagged / (double)newB.size()) >= Efficiency &&
        Max(Max(newB.extents(0),newB.extents(1)),newB.extents(2)) <= MaxWidth) {
      recurseOnThis.add(newB);
      return;
    } else {
      Cluster3R(flag,newB,Efficiency,MinWidth,MaxWidth,recurseOnThis);
      return;
    }
  }
  
  /* Find the slice positions */
  int istr=-1,   jstr=-1,   kstr=-1;
  int islice=-1, jslice=-1, kslice=-1;
  Cluster_Slice(i_sig,islice,istr,MinWidth*s[0]);
  Cluster_Slice(j_sig,jslice,jstr,MinWidth*s[1]);
  Cluster_Slice(k_sig,kslice,kstr,MinWidth*s[2]);

  // cout << "SLICE:"<<islice<<" "<<jslice<<" "<<kslice<<"\n";
  // cout << "STR  :"<<istr<<" "<<jstr<<" "<<kstr<<"\n";
  // cout.flush();

  /* Don't split small grids. */
  /* Bisectioning down to the minimum can occur, if a grid does not have any   */
  /* zero signature and the strength of its laplacian is lower than ITHRES. */
  /* We accept this grid, if its efficiency is greater than MIN_EFFICIENCY. */
  if ((islice+jslice+kslice==-3) ||
      (istr<MAXPOINTS2 && jstr<MAXPOINTS2 && kstr<MAXPOINTS2 &&
       (double)flagged/(double)flagbbox.size()>MIN_EFFICIENCY &&
       Max(Max(flagbbox.extents(0),flagbbox.extents(1)),flagbbox.extents(2)) 
       <= MaxWidth)) {
    // cout << "RETURNING " << flagbbox << endl << flush;
    recurseOnThis.add(flagbbox);
    return;
  }

  /* Split along the line with the biggest strength */
  BBox divme = flagbbox;
  if (istr > jstr && istr > kstr) {
    l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = divme.lower(2);
    u[0] = islice; u[1]= divme.upper(1); u[2] = divme.upper(2);
    tmpresult.add(BBox(3,l,u,s));
    l[0] = islice+s[0]; l[1]= divme.lower(1); l[2] = divme.lower(2);
    u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = divme.upper(2);
    tmpresult.add(BBox(3,l,u,s));
  } else if (jstr > kstr) {
    l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = divme.lower(2);
    u[0] = divme.upper(0); u[1]= jslice; u[2] = divme.upper(2);
    tmpresult.add(BBox(3,l,u,s));
    l[0] = divme.lower(0); l[1]= jslice+s[1]; l[2] = divme.lower(2);
    u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = divme.upper(2);
    tmpresult.add(BBox(3,l,u,s));
  } else {
    l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = divme.lower(2);
    u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = kslice;
    tmpresult.add(BBox(3,l,u,s));
    l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = kslice+s[2];
    u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = divme.upper(2);
    tmpresult.add(BBox(3,l,u,s));
  }

  /* foreach New bbox */
  for (BBox* newB = tmpresult.first(); newB; newB = tmpresult.next()) {
    /* Figure out how many points are flagged */
    // cout << "Handling "<< *newB <<"\n";
    flagged = 0;
    double myEff;
    for_3(ii, jj, kk, *newB, newB->stepsize())
      //cout << "SGF: "<<ii<<" "<<jj<<" "<<kk<<" "<<subgridflag(ii,jj,kk)<<"\n";
      //cout.flush();
      if (FastIndex3(flag,ii,jj,kk)>0) flagged++;
    end_for
    myEff = (double)flagged/(double)((*newB).size());
    /* Do we have ANY? */
    if (flagged) {
      // cout <<"Flagged " << flagged << " "<< *newB << "\n";
      // cout <<"Efficiency "<<myEff<<" of " << (*newB).size() <<"\n";
      // cout.flush();
      /* Stopping Condition */
      if ((myEff >= Efficiency && 
           Max(Max(newB->extents(0),newB->extents(1)),newB->extents(2)) <= MaxWidth) ||
          (newB->size() < MinWidth*MinWidth*MinWidth) ||
          (islice+jslice+kslice==-3)) {
        /* No. This one stays = the list if it contains flagged points */
        /* Don't add a degenerate grid! */
        int add = 1;
        for (BBox *tmp = recurseOnThis.first(); tmp; 
             tmp = recurseOnThis.next()) {
          if ((*tmp).inside((*newB).lower()) &&
              (*tmp).inside((*newB).upper())) add = 0;
        }
        if (add) {
          // cout << "  Adding "<< *newB <<"\n";
          // cout.flush();
          recurseOnThis.add(*newB);
        }
      } else {
        /* Yes: Recluster this sub-box */
        // cout <<"FALL: Flagged " << flagged << " "<< *newB << "\n";
        // cout <<"FALL: Efficiency "<<myEff<<" of " << (*newB).size() <<"\n";
        // cout.flush();
        BBoxList subcluster;
        Cluster3R(flag,*newB,Efficiency,MinWidth,MaxWidth,subcluster);
        for (BBox *me = subcluster.first(); me; me = subcluster.next()) 
          recurseOnThis.add(*me);
      }
    }
  }
  EndFastIndex1(i_sig);
  EndFastIndex1(j_sig);
  EndFastIndex1(k_sig);
  EndFastIndex3(flag);
  return;
}

void Cluster2(GridData(2)<short> &flag, 
              BBoxList& bblexclude, 
              double Efficiency,
              int MinWidth,
              int MaxWidth,
              int BufferWidth,
              BBoxList &Result) {
  BBoxList Result3D;
  GridData(2)<short> newflag(flag.bbox());
  BBox flat_3d(3,flag.bbox().lower(0),flag.bbox().lower(1),0,
               flag.bbox().upper(0),flag.bbox().upper(1),0,
               flag.bbox().stepsize(0),flag.bbox().stepsize(1),1);
  GridData(3)<short> clust_this(flat_3d);
  int flagged_points = 0;
  newflag = 0; clust_this = 0;

  for_2(i, j, flag.bbox(), flag.stepsize())
    if (flag(i,j)) {
      flagged_points ++;
      BBox where(2,i,j,i,j,
                 flag.bbox().stepsize(0),
                 flag.bbox().stepsize(1));
      where.grow(BufferWidth);
      newflag.equals(1,where);
    }
  end_for

  for (BBox *bbex=bblexclude.first();bbex;bbex=bblexclude.next())
    newflag.equals(0,*bbex);

  if (flagged_points) {
    for_2(i,j,flag.bbox(),flag.stepsize()) 
      clust_this(i,j,0) = newflag(i,j);
    end_for;
    Cluster3R(clust_this, clust_this.bbox(), Efficiency, 
              MinWidth, MaxWidth, Result3D);
  }

  /*
    Make sure that rank is appropriately set for the bboxes in 
    Result!!
    
    Manish Parashar Thu May  9 07:54:32 CDT 1996
    */
  for (BBox* tmpbb=Result3D.first();tmpbb;tmpbb=Result3D.next()) {
    BBox tmp2d(2,tmpbb->lower(),tmpbb->upper(), tmpbb->stepsize());
    Result.add(tmp2d);
  }
}

void Cluster3(GridData(3)<short> &flag, 
              BBoxList& bblexclude, 
              double Efficiency,
              int MinWidth,
              int MaxWidth,
              int BufferWidth,
              BBoxList &Result) {
  BBoxList tmpresult, result;
  GridData(3)<short> newflag(flag.bbox());
  int flagged_points = 0;
  newflag = 0;
  
  for_3(i, j, k, flag.bbox(), flag.stepsize())
    if (flag(i,j,k)) {
      flagged_points ++;
      BBox where(3,i,j,k,i,j,k,
                 flag.bbox().stepsize(0),
                 flag.bbox().stepsize(1),
                 flag.bbox().stepsize(2));
      where.grow(BufferWidth);
      newflag.equals(1,where);
    }
  end_for

  for (BBox *bbex=bblexclude.first();bbex;bbex=bblexclude.next()) 
    newflag.equals(0,*bbex);
  
  if (flagged_points) {
    Cluster3R(newflag, newflag.bbox(), Efficiency, 
              MinWidth, MaxWidth, Result);
  }
}

---

Quickstart     Users Guide     Programmers Reference     Installation      Examples     Download

---

AMROC Main      Home      Contact

last update: 06/01/04