How to use modified Shepard method

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

I've looked at your code, there is an error in the interpolation part.

You wrote:

Code:

//////////////////////////////////////////////////////
// calc interpolant
//////////////////////////////////////////////////////
for( int i = 0; i < n; i++ )
{
    x(0) = xy(i, 2);
    //ap::vmove(&x(0), 1, &xy(i, 2), 1, ap::vlen(0,nx-1));
    double val = idwcalc(z1, x);
    idwerrors = ap::fp_neq(val, xy(i,nx));
}

But x was not correctly initialized. Correct code will be

Code:

//////////////////////////////////////////////////////
// calc interpolant
//////////////////////////////////////////////////////
for( int i = 0; i < n; i++ )
{
    x(0) = xy(i, 0);
    x(1) = xy(i, 1);
     double val = idwcalc(z1, x);
    idwerrors = ap::fp_neq(val, xy(i,nx));
}

So interpolant will match function value at the points. However, you should also look at the interpolation quality between points, because IDW sometimes gets tricky.

AlekseyL · **Joined:** Mon Sep 13, 2010 12:58 am **Posts:** 14

Hi Sergey. Thanks for reply. ;)

I fixed my code. And now idwcalc returns the same values which I initialize my xy array. So I have question - how can I get approximation values instead zero values in my data (see signal.txt in IDWInterpolation.zip)?

Best regards, Aleksey.

P.S. I'm start to doubt that I chose right way for my task. :(

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Quote:

And now idwcalc returns the same values which I initialize my xy array. So I have question - how can I get approximation values instead zero values in my data (see signal.txt in IDWInterpolation.zip)?

Well, you've specified zero values at these points, so it returns zero values. If you don't know function value at these points, simply don't pass them to interpolant construction function. You can pass less then width*height points, and they don't have to be on a regular grid. I also recommend to use IDWBuildModifiedShepardR() function.

However, I doubt that you will get adequate results with any kind of interpolation in your case. You have very thin strip of points with known values and want to extend it to the whole plane. But there is too little information to reconstruct function surface far from the diagonal of your matrix. You can expect good results near the diagonal, between diagonal points, but not far away from them.

AlekseyL · **Joined:** Mon Sep 13, 2010 12:58 am **Posts:** 14

Thanks for reply!

I don't have good news again :(
I've changed my code as you say. I've inited xy array only with non zero values:

Code:

int k = 0;
   for ( int i = 0; i < nArrayHeight; i++ )
   {
      for ( int j = 0; j < nArrayWidth; j++ )
      {
         xy(k, 0) = i;
         xy(k, 1) = j;
         if ( m_ppArray[i][j].dValue != 0 )
         {
              xy(k, 2) = m_ppArray[i][j].dValue;
         }
         k++;
      }
   }

And as you recommended I've used IDWBuildModifiedShepardR() function:

Code:

int r = 3;

   //idwbuildmodifiedshepard(xy, n, nx, d, nq, nw, z1);
   idwbuildmodifiedshepardr(xy, n, nx, r, z1);

But now idwcalc returns the same values as inited points and all another values are equal to -6.2774385622041925e+066:

Code:

for( int i = 0; i < n; i++ )
   {
      x(0) = xy(i, 2);
      double val = idwcalc(z1, x);
   }

Where I did mistake again?

Thanks
Best regards, Aleksey.

P.S. In attached file whole code of method.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

You should write it as follows:

Code:

int k = 0;
for ( int i = 0; i < nArrayHeight; i++ )
{
    for ( int j = 0; j < nArrayWidth; j++ )
        if ( m_ppArray[i][j].dValue != 0 )
        {
            xy(k, 0) = i;
            xy(k, 1) = j;
            xy(k, 2) = m_ppArray[i][j].dValue;
            k++;
        }
}

and pass k as the number of points.

Your original code initializes only points with nonzero function values, but passes random mix of initialized and uninitialized points.

AlekseyL · **Joined:** Mon Sep 13, 2010 12:58 am **Posts:** 14

Now I initialize my xy array as you say. Then call method

Code:

idwbuildmodifiedshepardr(xy, k, nx, r, z1);

with k points.

Maybe should I change this code? :

Code:

x.setlength(2);
for( int i = 0; i < n; i++ )
{
   x(0) = xy(i, 0);
   x(1) = xy(i, 1);
   double val = idwcalc(z1, x);
}

Because now idwcalc returns first k values the same which was initialized in xy and all other equal to 74134.199999999997.

Thanks
Best regards, Aleksey.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Quote:

Because now idwcalc returns first k values the same which was initialized in xy and all other equal to 74134.199999999997.

74134 is returned for off-diagonal points, right? It should be related to the inability of interpolant to extrapolate function beyond known set. As I told you, your problem seems insolvable - you have only thin strip of points with known values and want to extrapolate it too far away. You just don't have enough data to do such extrapolation.

Of course, you can have some problem-specific assumptions about properties of your function, but in such case you have to write problem-specific interpolation algorithm.

AlekseyL · **Joined:** Mon Sep 13, 2010 12:58 am **Posts:** 14

Hi Sergey.

Yes you are right, for all zero value points idwcalc returns 74134.

But if you'll watch to my data set (see attached file please), you can see that there is another one "diagonal ", so I think that it able to calc interpolation values between "diagonals". And I don't need calc all zero value points, even though 2-3 near points for my known points.

So if you certain, that alglib IDW methods will not help me, could you recommend me interpolation algorithm which could lick my problem, please?

Thanks for your time
Best regards, Aleksey.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

You can try increasing R, influence radius. It must be large enough to make any point between two diagonals influenced by both diagonals. I think, that R=10 or R=15 will be enough.

As for other methods, you can look at the thin plate splines or penalized bicubic splines. However, ALGLIB doesn't includes any of these methods yet.

AlekseyL · **Joined:** Mon Sep 13, 2010 12:58 am **Posts:** 14

I've tried to increase R value to 15 and then 20, but the result is the same.

Ok Sergey, thank you for your spending time. I hope this topic will useful for another one who want to use Shepard method from your lib.

All the best! ;)
Best regards, Aleksey.

forum.alglib.net

Forum rules

How to use modified Shepard method

Who is online