2D '& 3D Surface Contour Map



Click here for a larger image.



Click here for a larger image.

Environment: VC6

A year and a half year ago, I published this article to the Codeguru site and got a number of requests about the Kriging algorithm contour map. Unfortunately, my project was changed shortly after that article and later I quit the company so I couldn't find time to finish this Contour business. A week ago, I happened to need a contour map again so I decided to solve the Kriging algorithm. I searched the Internet for a commercial library but they all look ugly and hard to use. So, I made up my mind to make my own algorithm. The Kriging algorithm is easy to find, but this algorithm needs a Matrix and solver (LU-Decomposition). Again, I couldn't find suitable code for this. I tried to use GSL first but this made my code too big and was slower. Finally, I went back to "Numerical Recipe in C"—yes, that horrible-looking C code—and changed the code there to my taste.

If you read this article before, the rendering part hasn't been changed much. I added the Kriging algorithm and revised the codes a little bit. Following is the Kriging Algorithm:

template<class ForwardIterator>
double GetDistance(const ForwardIterator start, int i, int j)
{
  return ::sqrt(::pow(((*(start+i)).x - (*(start+j)).x), 2) +
                ::pow(((*(start+i)).y - (*(start+j)).y), 2));
}

template<class ForwardIterator>
double GetDistance(double xpos, double ypos,
                   const ForwardIterator start, int i)
{
  return ::sqrt(::pow(((*(start+i)).x - xpos), 2) +
                ::pow(((*(start+i)).y - ypos), 2));
}

template<class T, class ForwardIterator>
class TKriging : public TInterpolater<ForwardIterator>
{
public:
  TKriging(const ForwardIterator first, const ForwardIterator last,
           double dSemivariance) : m_dSemivariance(dSemivariance)
  {
    m_nSize = 0;
    ForwardIterator start = first;
    while(start != last) {
      ++m_nSize;
      ++start;
    }

    m_matA.SetDimension(m_nSize, m_nSize);

    for(int j=0; j<m_nSize; j++) {
      for(int i=0; i<m_nSize; i++) {
        if(i == m_nSize-1 || j == m_nSize-1) {
          m_matA(i, j) = 1;
          if(i == m_nSize-1 && j == m_nSize-1)
            m_matA(i, j) = 0;
          continue;
        }
        m_matA(i, j) = ::GetDistance(first, i, j) * dSemivariance;
      }
    }
    int nD;
    LUDecompose(m_matA, m_Permutation, nD);
  }
  double GetInterpolatedZ(double xpos, double ypos,
                          ForwardIterator first,
                          ForwardIterator last)
                          throw(InterpolaterException)
  {
    std::vector<T> vecB(m_nSize);
    for(int i=0; i<m_nSize; i++) {
      double dist = ::GetDistance(xpos, ypos, first, i);
      vecB[i] = dist * m_dSemivariance;
    }
    vecB[m_nSize-1] = 1;

    LUBackSub(m_matA, m_Permutation, vecB);

    double z = 0;
    for(i=0; i<m_nSize-1; i++) {
      double inputz = (*(first+i)).z;
      z += vecB[i] * inputz;
    }
    if(z < 0)
      z = 0;
    return z;
  }
private:
  TMatrix<T> m_matA;
  vector<int> m_Permutation;
  int m_nSize;
  double m_dSemivariance;
};

typedef TKriging<double, Point3D*> Kriging;

Because of the template, this doesn't look that clean but you can get the idea if you look at it carefully. The matrix solver is as follows:

template<class T>
void LUDecompose(TMatrix<T>& A, std::vector<int>&
                 Permutation, int& d) throw(NumericException)
{
  int n = A.GetHeight();
  vector<T> vv(n);
  Permutation.resize(n);

  d=1;

  T amax;
  for(int i=0; i<n; i++) {
    amax = 0.0;
    for(int j=0; j<n; j++)
      if(fabs(A(i, j)) > amax)
        amax = fabs(A(i, j));

    if(amax < TINY_VALUE)
      throw NumericException();

    vv[i] = 1.0 / amax;
  }

  T sum, dum;
  int imax;
  for(int j=0; j<n; j++) {
    for (i=0; i<j; i++) {
      sum = A(i, j);
      for(int k=0; k<i; k++)
        sum -= A(i, k) * A(k, j);
      A(i, j) = sum;
    }
    amax = 0.0;

    for(i=j; i<n; i++) {
      sum = A(i, j);
      for(int k=0; k<j; k++)
        sum -= A(i, k) * A(k, j);

      A(i, j) = sum;
      dum = vv[i] * fabs(sum);

      if(dum >= amax) {
        imax = i;
        amax = dum;
      }
    }

    if(j != imax) {
      for(int k=0; k<n; k++) {
        dum = A(imax, k);
        A(imax, k) = A(j, k);
        A(j, k) = dum;
      }
      d = -d;
      vv[imax] = vv[j];
    }
    Permutation[j] = imax;

    if(fabs(A(j, j)) < TINY_VALUE)
      A(j, j) = TINY_VALUE;

    if(j != n) {
      dum = 1.0 / A(j, j);
      for(i=j+1; i<n; i++)
        A(i, j) *= dum;
    }
  }
}

template<class T>
void LUBackSub(TMatrix<T>& A, std::vector<int>&
               Permutation, std::vector<T>& B)
               throw(NumericException)
{
  int n = A.GetHeight();
  T sum;
  int ii = 0;
  int ll;
  for(int i=0; i<n; i++) {
    ll = Permutation[i];
    sum = B[ll];
    B[ll] = B[i];
    if(ii != 0)
      for(int j=ii; j<i; j++)
        sum -= A(i, j) * B[j];
    else if(sum != 0.0)\
      ii = i;
    B[i] = sum;
  }

  for(i=n-1; i>=0; i--) {
    sum = B[i];
    if(i< n) {
      for(int j=i+1; j<n; j++)
        sum -= A(i, j) * B[j];
    }
    B[i] = sum / A(i, i);
  }
}

By using this algorithm, making a 3D grid is easy. Let's assume we're making a 200x200 grid and we have some scattered data. Then, what we need to do is this:

  vector<Point3D> input    // assume this vector has KNOWN 3D points

  Interpolater* pInterpolater = new Kriging(input.begin(),
                                    input.end(), 4);

  vector<double> vecZs;

  for(int j=0; j<200; j++) {
    for(int i=0; i<200; i++) {
      vecZs.push_back(pInterpolater->GetInterpolatedZ(i, j,
                                                      input.begin(),
                                                      input.end()));
    }
  }
  // Now, vecZs has 40000 z values

  delete pInterpolater;

If you have all the grid points with 3D data, you can make a bitmap file with it, or make a triangle strip to render with OpenGL. If you remember that the old contour map was produced from an InverseDistanced algorithm (you can switch to Inverse Distance in the Option menu), you'll find a vast improvement over it. I compared the Kriging generated contour map with some commercial programs, and they were almost identical. I hope this helps programmers who want to make a contour map.

Downloads

Download demo project - 221 Kb
To use this project, open one of the "*.txt" files first, then click the "2D" or "3D" button in the toolbar.



Comments

  • together ancestor lay across boyfriend number physical

    Posted by kjcgzyqk on 06/13/2013 11:11am

    nationality Longchamp Bags Oultet remember goldfish simple water stormy franc school modal float twentieth Longchamp UK glad willing sleeve importance newspaper Instyler Hair Straightener ocean bottom loose grassland laugh bear hire gathering mine imperialist boyfriend http://www.insanityuk.net/ user company frog midday

    Reply
  • isabel marant sneakers isabelsneakersstore.slicins.com

    Posted by INSUNDELENT on 06/08/2013 06:45pm

    14: 45-89Choosing Health A consultation on action to improve peoples health, Spring 2004. isabel marant

    Reply
  • isabel marant sneakers

    Posted by demGredybom on 03/07/2013 02:49pm

    isabel marant sneakers trail select wheelchair rubbish mean dumpling safe pie trainer place isabel marant sneakers round warmhearted can awake Britain expedition how travel object midnight isabel marant shoes devotion conduct schoolgirl non-identical generally smooth chime soap everyday verb Isabel Marant Dicker Suede Ankle Boots split anger research object major gentleman carry exercise satisfaction national isabel marant ankle boots immediately listener worthy gather arisen coral work Russian slavery fast isabel marant sneakers sale container important talk October robot nurse both lightly poster cloth http://isabelmarantstore0035.webs.com balloon politician dull http://isabelmarantstore0040.webs.com get most

    Reply
  • nike free run 2

    Posted by diomormFelo on 03/03/2013 02:36pm

    nike free run main exam own delay fully invitation merciful faith jump that nike free run 2 sow Atlantic unless job unreserved thunder obey assistant best-seller defence nike free run 3 air packet farmer instruction diamond garden happen sorry personally gray nike free run black handtruck fuel class receive safe ordinary senior thick ball until nike free sale Wednesday exact carpet color hair thirteen nation care troop inspiring nike free run row extra through dislike strange mountain India laziness snow locust http://nikefreerun124.webs.com waiter situation funny http://nikefreerun124.webs.com each researcher

    Reply
  • nike free run 3

    Posted by gegeneume on 03/03/2013 11:16am

    nike free run purpose expression pupil clever section come hilly spit institute best-seller nike free run 2 pie speed mad salty fiercly judge place button argument bad nike free run 3 message wherever handkerchief swift pottery sixteen particular world advertisement fitting nike free run black clear child oh hillside tank whisper these twin sing clerk nike free shoes packet village window basket summary weight architechure abroad suffer green nike free run living-room finally remain interview energy plough open tape-recording quiet brotherhood http://nikefreerun128.webs.com alike spaceship win http://nikefreerun124.webs.com magazine determination

    Reply
  • nike free run womens

    Posted by TredHeava on 03/02/2013 04:38pm

    nike free run ring job Japanese key thermos shallow coral thirtieth exact granny nike free run 2 dye explode cleaning burst pence tax lamb side sadly centimetre nike free run 3 hawk high-rise wire feeling exploit moreover sweater make history short nike free run womens custom storage experiment she Swiss review rockface hatch pack policeman nike free shoes sure party describe gain his outside get beam meat lack nike free run zone button toilet decision moneylender mourn difference biscuit relativity wall http://nikefreerun126.webs.com toward industry polite http://nikefreerun1233.webs.com flashlight chocolate

    Reply
Leave a Comment
  • Your email address will not be published. All fields are required.

Top White Papers and Webcasts

  • Live Event Date: May 6, 2014 @ 1:00 p.m. ET / 10:00 a.m. PT While you likely have very good reasons for remaining on WinXP after end of support -- an estimated 20-30% of worldwide devices still are -- the bottom line is your security risk is now significant. In the absence of security patches, attackers will certainly turn their attention to this new opportunity. Join Lumension Vice President Paul Zimski in this one-hour webcast to discuss risk and, more importantly, 5 pragmatic risk mitigation techniques …

  • The impact of a data loss event can be significant. Real-time data is essential to remaining competitive. Many companies can no longer afford to rely on a truck arriving each day to take backup tapes offsite. For most companies, a cloud backup and recovery solution will eliminate, or significantly reduce, IT resources related to the mundane task of backup and allow your resources to be redeployed to more strategic projects. The cloud - can now be comfortable for you – with 100% recovery from anywhere all …

Most Popular Programming Stories

More for Developers

Latest Developer Headlines

RSS Feeds