Shoeisha Technology Media

CodeZine(コードジン)

記事種別から探す

離散した点を補間してグラフを描画する方法

Lagrange、Newton、Splineなどの各種補間方式の実装

  • LINEで送る
  • このエントリーをはてなブックマークに追加
2005/09/16 12:00

実験結果をグラフに表示する時、測定点をプロットしてから、その間を自在定規などで結んだ経験があると思います。本稿では、これを自動的に行うプログラムを紹介します。このような作業は「補間」(Interpolation)と呼ばれ、色々な方式があります。

はじめに

 実験結果をグラフに表示する時、測定点をプロットしてから、その間を自在定規などで結んだ経験があると思います。ここでは、これを自動的に行うプログラムを紹介します。このような作業は「補間」(Interpolation)と呼ばれますが、「補間」には色々な方式がありますので、これらを比較できるようにしました。

対象読者

 LagrangeNewtonSplineなどの補間方式に興味を持ち、実験結果をまとめるのに、自分で作ったプログラムを応用したい人。また、CGの基礎としての補間の原理を学びたい人。

必要な環境

 J2SE 5.0を使っていますが、それより古いバージョンでも大丈夫です。

補間とは

 「狭義の補間」は、二点のデータを正確なものと仮定し、その中間点のデータを推測することで、「内挿」とも言います。これは、デジカメなどの画像処理で用いられます。一方、「広義の補間」は、与えられたすべてのデータ点を通るように曲線を描くこと、および、その曲線の方程式を求めることで、データ点に誤差があっても、そのまま補間されますので、補間は正しい値を推測するものではありません。補間は、実験結果のグラフ作成のほか、アウトラインフォントCGなどに広く用いられています。

 似たものに、「曲線の当てはめ(Curve Fitting)」があります。これは、あらかじめ与えられた次数の多項式で、測定点との誤差を最小にするもので、実験結果のグラフ化や実験式の作成に用いられます。この場合は、測定点を必ず通るのではなく、より正しい値を推測します。

補間のいろいろと、プログラムでの実現方法

 0番目からn番目まで(n+1)個の点を必ず通る曲線を表現するには、n次の多項式が必要十分条件で、この多項式を求めることが補間です。

 多項式を求める方法には、Lagrange多項式が良く知られていますが、Newtonの方法もあり、二つの方法の結果は同じになります。後者は、データ点の追加が容易な長所がある代わりに、再帰処理のために時間が掛かります。筆者の案によるLagrange補間の改良(すでに誰かがやっているかも知れませんが)も、併せて紹介します。

 一方、多項式の次数を限定し、狭い区間ごとに異なる多項式を求める方法にSplineがあります。これは、曲線の円滑な連続性には欠けますが、異常な振動が起きません。Splineの中では、ここで紹介する3次Splineが代表的なものです。CGなどで知られるNURBSは、Non Uniform Rational B-Splinesの略で、高度なSpline補間方式が用いられています。

線形補間

 相隣る二点間を直線(一次式)で結ぶ一番単純な補間で、学校のレポートでグラフを描く時には、「最もいけない」とされるものです。プログラム的にも面白味がありませんので、ここでは取り上げません。

Lagrange補間

 データ点を(X0,Y0)、(X1,Y1)、…、(Xn,Yn)の(n+1)個とすると、これらの点を必ず通るn次の多項式Pn(x)は下式で求まります。LはLagrange、PはPolynomial(多項式)の意味です。

 プログラムでは、データの個数をdata_countとしていますので、j0からdata_count-1の間で変化させ、

(x-dataX[j])/(dataX[i]-dataX[j])

 を掛け合わせて結果をi番目のlagとし、ylag*dataY[i]を加算して求めます。

Lagrange補間の改良

 上記のLagrange補間は、両端付近で大きな異常振動が生じる場合があります。これを緩和するために、次の改良を行います。

  1. データ点の(X0,Y0)、(X1,Y1)、…、(Xn,Yn)の(n+1)個に対して、両端に合計2個の点を追加します。入力した(X0,Y0)、(X1,Y1)、…、(Xn,Yn)を、(X1,Y1)、(X2,Y2)、…、(Xn,Yn)、(Xn+1,Xn+1)のように名前を付け替え、両端に新たに(X0,Y0)と(Xn+2,Yn+2)を追加します。プログラムでは、名前の付け替えは、
  2. dataX1[i+1]=dataX[i] と dataY1[i+1]=dataY[i] 
    
    によります。
  3. (X0,Y0)は、(X1,Y1)と(X2,Y2)を線形的に外挿して求め、(Xn+2,Yn+2)は、(Xn,Yn)と(Xn+1,Yn+1)を線形的に外挿します。
  4. プログラムでは、
    dataX1[0]=2*dataX[0]-dataX[1]
    dataX1[data_count+1]=2*dataX[data_count-1]-dataX[data_count-2]
    
    などによります。
  5. 以上の(n+3)個の点に対して通常のLagrange補間を行いますが、追加した(X0,Y0)と(X1,Y1)の間、および(Xn,Yn)と(Xn+1,Yn+1)の間は表示しません。
  6. プログラムでは、
    if(x>=dataX1[1] && x<=dataX1[data_count])
    
    の条件を付けて描画します。

 実行結果を見ると、ある程度の効果があることが分りますが、線形的な外挿が必ずしも好ましくない場合があり、かえって振動を増加させることがあります。

Newton補間

 データ点を(X0,Y0)、(X1,Y1)、…、(Xn,Yn)の(n+1)個とすると、これらの点を必ず通るn次の多項式Pn(x)は下式で得られます。Ciは差分商(Divided Difference) と呼ばれ、〔XiXi-1....X0〕で表わしますが、この計算は、再帰プログラムで容易に実現できます。

 プログラムでは、n[i](x-dataX[j])の乗算(ただし、j0からi-1の間)で求めます。

 c[i]は、メソッドcalculateDividedDiff(i,0)で求め、yは、c[i]*n[i]の加算で求めます。

Spline補間

 Spline補間は、末端部の条件や区間の取り方により何種類かありますが、結果はほとんど同じになります。データ点(Xi-1,Yi-1)と(Xi,Yi)の間を結ぶSpline補間の一例であり、広く用いられている3次Splineの式を下記に示します。

 ただし、

 プログラムでは、3次Splineの第1項から第4項までを、yy0yy1yy2yy3として計算し、加算します。

 一次微分はdif1[i]、二次微分はdif2[i]で表わし、上式をそのままコード化します。

プログラム

import java.applet.Applet;
import java.awt.*;
import java.awt.event.*;

public class Interpolation extends Applet implements 
   MouseListener,MouseMotionListener,ActionListener{

   static final int X0=50,Y0=50;   //「方眼紙」の左上座標
   int data_count;                 //データ個数

   //入力データのx値(座標値ではない)
   double[] dataX=new double[20];
   //入力データのy値(座標値ではない)
   double[] dataY=new double[20]; 

   Button[] bt=new Button[6];

   public void init(){

      data_count=0;

      // ---------------------- マウス関係準備 ----------------------

      addMouseMotionListener(this); 
      addMouseListener(this);

      // ---------------------- ボタン関係準備 ----------------------

      String[] label={"  Lagrange  ","   Newton   ",
                      " Lagrange改 "," 3次Spline ",
                      "  補間消去  "," データ消去 "};

      for(int i=0;i<6;i++){
         bt[i]=new Button(label[i]);
         add(bt[i]);
         bt[i].addActionListener(this);
      }

      //ボタンの色を設定する
      bt[0].setBackground(Color.green);   //Lagrange
      bt[1].setBackground(Color.yellow);  //Newton
      bt[2].setBackground(Color.magenta); //Lagrange改
      bt[3].setBackground(Color.cyan);    //3次Spline

   }

   // ------------------- マウス移動関係メソッド --------------------

   public void mouseMoved(MouseEvent me){  //マウスを動かした

      int x=me.getX();
      int y=me.getY();

      float xx=(float)((x-X0)/10.0);   //floatにキャストして丸める
      float yy=(float)(40.0-(y-Y0)/10.0);

      Graphics g=getGraphics();

      g.setColor(Color.white);
      g.fillRect(X0-5,Y0+445,300,20);  //前の表示を消す

      g.setColor(Color.black);
      if(xx>=0.0 && xx<=60.0 && yy>=0.0 && yy<=40.0)
         // 座標の表示
         g.drawString("X座標= "+ xx +" ,Y座標= "+ yy ,X0,Y0+460);
      else
         g.drawString("範囲外です。",X0,Y0+460);

   }

   //使用しないが、記述しておかないとエラーになる
   public void mouseDragged(MouseEvent me){ }

   // ------------------- マウス動作関係メソッド --------------------

   public void mousePressed(MouseEvent me){  //マウスをクリックした

      int x=me.getX();
      int y=me.getY();

      double xx=(x-X0)/10.0;
      double yy=40.0-(y-Y0)/10.0;

      if(xx>=0.0 && xx<=60.0 && yy>=0.0 && yy<=40.0){

         dataX[data_count]=xx;
         dataY[data_count]=yy;
         data_count++;
         repaint();

      }

   }

   //使用しないが、記述しておかないとエラーになる
   public void mouseClicked(MouseEvent me){ }
   public void mouseEntered(MouseEvent me){ }
   public void mouseReleased(MouseEvent me){ }
   public void mouseExited(MouseEvent me){ }

   // --------------------- ボタン関係メソッド ----------------------

   //ボタンがクリックされた
   public void actionPerformed(ActionEvent ae){

      if(ae.getSource()==bt[0])    drawLagrange();    //Lagrange
      if(ae.getSource()==bt[1])    drawNewton();      //Newton
      if(ae.getSource()==bt[2])    drawLagrange1();   //Lagrange改
      if(ae.getSource()==bt[3])    drawSpline();      //3次Spline
      if(ae.getSource()==bt[4])    repaint();         //補間消去
      if(ae.getSource()==bt[5])    clearData();       //データ消去 

   }

   // --------------------- 補間処理メソッド類 ----------------------

   //Lagrange補間
   private void drawLagrange(){
      
      double x,y,lag;

      Graphics g=getGraphics();
      g.setColor(new Color(0,192,0));   //dark greenで線を引く

      for(x=dataX[0];x<=dataX[data_count-1];x+=0.01){
         y=0.0;
         for(int i=0;i<data_count;i++){
            lag=1.0;
            for(int j=0;j<data_count;j++)
               if(j!=i) lag*=(x-dataX[j])/(dataX[i]-dataX[j]);
            y+=lag*dataY[i];
         }
         
         g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40.0-y)),0,0);
 
      }

   }

   //Newton補間
   private void drawNewton(){

      int i,j;
      double x,y,p;
      double[] n=new double[20];
      double[] c=new double[20];

      Graphics g=getGraphics();
      g.setColor(new Color(192,192,0));   //dark yellowで線を引く

      for(x=dataX[0];x<=dataX[data_count-1];x+=0.01){
         n[0]=1.0;
         for(i=1;i<data_count;i++){
            p=1.0;
            for(j=0;j<i;j++)
               p*=x-dataX[j];
            n[i]=p;
         }

         for(i=0;i<data_count;i++)
            c[i]=calculateDividedDiff(i,0);

         y=0.0;
         for(i=0;i<data_count;i++)
            y+=c[i]*n[i];

         g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40.0-y)),0,0);

      }

   }

   //Newton補間から呼び出される再帰メソッド
   private double calculateDividedDiff(int a,int b){

      if(a<=b) return dataY[a];
      else     return (calculateDividedDiff(a,b+1)-
         calculateDividedDiff(a-1,b))/(dataX[a]-dataX[b]);

   }

   //Lagrange補間の改良
   private void drawLagrange1(){

      double x,y,lag;
      double[] dataX1=new double[22];    //データ点を増やす準備
      double[] dataY1=new double[22];

      for(int i=0;i<data_count;i++){     //入力したデータ点を移す
         dataX1[i+1]=dataX[i];
         dataY1[i+1]=dataY[i];
      }

      //両端のデータを外挿して求める
      dataX1[0]=2*dataX[0]-dataX[1];
      dataY1[0]=2*dataY[0]-dataY[1];
      dataX1[data_count+1]=2*dataX[data_count-1]-dataX[data_count-2];
      dataY1[data_count+1]=2*dataY[data_count-1]-dataY[data_count-2];
      
      Graphics g=getGraphics();
      g.setColor(new Color(192,0,192));   //dark magentaで線を引く

      for(x=dataX1[0];x<=dataX1[data_count+1];x+=0.01){
         y=0.0;
         for(int i=0;i<data_count+1;i++){
            lag=1.0;
            for(int j=0;j<data_count+1;j++)
               if(j!=i) lag*=((x-dataX1[j])/(dataX1[i]-dataX1[j]));
            y+=(dataY1[i]*lag);
         }

         //実際に入力した範囲のみを描画
         if(x>=dataX1[1] && x<=dataX1[data_count])
            g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40.0-y)),0,0);
 
      }

   }

   //3次Spline補間
   private void drawSpline(){

      int i;
      double x,y,yy0,yy1,yy2,yy3;
      double[] h=new double[20];     //間隔
      double[] dif1=new double[20];  //一次微分
      double[] dif2=new double[20];  //二次微分

      h[0]=0.0;
      dif2[0]=0.0;
      dif2[data_count-1]=0.0; 

      Graphics g=getGraphics();
      g.setColor(new Color(0,192,192));  //dark cyanで線を引く

      for(i=1;i<data_count;i++){
         h[i]=dataX[i]-dataX[i-1];             //間隔を計算
         dif1[i]=(dataY[i]-dataY[i-1])/h[i];   //一次微分を計算
      }

      for(i=1;i<data_count;i++)
         //二次微分を計算
         dif2[i]=(dif1[i+1]-dif1[i])/(dataX[i+1]-dataX[i-1]);

      i=1;
      for(x=dataX[0];x<dataX[data_count-1];x+=0.01){
         if(x<dataX[i]){
            yy0=dif2[i-1]/(6*h[i])*(dataX[i]-x)
               *(dataX[i]-x)*(dataX[i]-x);       //第1項
            yy1=dif2[i]/(6*h[i])*(x-dataX[i-1])
               *(x-dataX[i-1])*(x-dataX[i-1]);   //第2項
            yy2=(dataY[i-1]/h[i]-h[i]*dif2[i-1]/6)
               *(dataX[i]-x);                    //第3項
            yy3=(dataY[i]/h[i]-h[i]*dif2[i]/6)*
               (x-dataX[i-1]);                   //第4項
            y=yy0+yy1+yy2+yy3;
            g.drawRect(X0+(int)(10*x),Y0+(int)(10*(40-y)),0,0);
         }
         else i++;
      }

   }

   //データ点消去
   private void clearData(){
      
      data_count=0;
      for(int i=0;i<20;i++){
         dataX[i]=0.0;
         dataY[i]=0.0;
      }
      repaint();

   }

   public void paint(Graphics g){

      int i,j,x,y;

      //「方眼紙」を描く
      g.setColor(Color.lightGray);
      for(i=0;i<=60;i++)
         g.drawLine(X0+10*i,Y0,X0+10*i,Y0+400);
      for(j=0;j<=40;j++)
         g.drawLine(X0,Y0+10*j,X0+600,Y0+10*j);
      g.setColor(Color.gray);
      for(i=0;i<=12;i++)
         g.drawLine(X0+50*i,Y0,X0+50*i,Y0+400);
      for(j=0;j<=8;j++)
         g.drawLine(X0,Y0+50*j,X0+600,Y0+50*j);

      //目盛を入れる
      g.setColor(Color.black);
      for(i=0;i<=12;i++)
         g.drawString(String.valueOf(5*i),X0-5+50*i,Y0+424);
      for(j=0;j<=8;j++)
         g.drawString(String.valueOf(5*(8-j)),X0-25,Y0+5+50*j);
      

      //データ点を描く
      g.setColor(Color.red);

      i=data_count-1;
      while(i>=0){
         x=X0+(int)(dataX[i]*10);
         y=Y0+(int)((40.0-dataY[i])*10);
         g.drawOval(x-2, y-2, 4, 4); 
         i--;
      }
       
   }

}

プログラムの使い方

 アプレットを実行すると、上部に、左から[Lagrange][Newton][Lagrange改][3次Spline][補間消去][データ消去]のボタンが表示され、下に方眼紙状のデータ表示領域が現われます。

 「方眼紙」内でクリックすると、データ点が入力され、赤丸が表示されます。マウスポインタの位置は、左下に数字で表わされますので、参考にします。データ入力は、必ず左から行ってください。途中で引き返すと、その点までしか表示されません(これは改善できますが、プログラムを簡単化するために省略しました)。間違ってデータを入力した場合は、[データ消去]のボタンをクリックして、最初からやり直してください。

 また、入力できるデータ点の数は20個に制限してありますので、注意してください。

 [Lagrange]などのボタンをクリックすると、補間が実行され、それぞれ、ボタンと同じ色で曲線が描かれます。次々とボタンをクリックすると、対応する補間曲線が追加されます。[補間消去]ボタンをクリックすると、データ点を残して補間曲線のみが消えるので、何回でも繰り返すことができます。

 データを入力し直すときは、[データ消去]をクリックしてください。

 なお、[Lagrange]と[Newton]は全く同じ結果になります。両者を実行できるようにしたのは、その事実の確認のためです。

プログラムの実行結果

 下図において、赤丸は、「暴れ」が出やすいように、やや意地悪に入力したデータ点(普通の実験データの場合は、きれいな曲線になります)、緑の線はLagrange、マゼンタはLagrange改良型、シアンは3次Splineによる補間結果です。Newton補間はLagrange補間と結果が同じなので、表示していません。Lagrangeに比べ、Lagrange改良型は暴れ方が減っていることがわかります。3次Splineは、「暴れ」はありませんが、線がスムーズでない部分があります。

 「方眼紙」の下は、マウスポインタが指している座標で、この例では、一番右のデータ点を指しています。

図 実行結果の一例
図 実行結果の一例

まとめ

 ここで紹介した補間プログラムは、プログラム的には比較的単純なものです。各種の補間方式の数学的意味が理解できれば、後は、それをコーディングするだけで、特別なアルゴリズムも要りません。プログラムが得意な人は、数学も得意でしょうから(筆者は少なくともそれを期待しています)、数学的な詳しい説明は省きました。

 このプログラムは、プログラムを作成する楽しみもありますが、実際にデータを入れてみて、どのような結果になるかを試してみる面白さがあります。電気に詳しい人であれば、補間結果の「暴れ」の原因がお分りかと思います。データ群は高い周波数成分を有するのに、強制的にn次式でしか表わされず、高調波成分がカットされるからです。

参考資料

 この記事は、筆者のウエブサイト『Visual C++の易しい使い方-Lagrange、Newton、Splineなどの色々な補間方式を実現する-(ただし、現在はVisual C++ 2005 Express Edition版に改定)』をJava言語に書き直し、分かりやすく加筆したものです。

 補間は、数値計算の基本で、初期のコンピュータの技術計算への応用の主なものの一つでした。したがって、

 などの数学辞典類にはたいてい載っています。プログラムに限れば、

  1. 『C & FORTRANによる数値解析の基礎』 川崎晴久 著、共立出版、1993年5月
  2. 『CとJavaで学ぶ数値シミュレーション入門』 峯村吉泰 著、森北出版、1999年4月
  3. 『ANSI Cによる数値計算法入門』 堀之内総一 他著、森北出版、2002年12月

 などがあります。



  • LINEで送る
  • このエントリーをはてなブックマークに追加

修正履歴

  • 2008/03/16 10:55 参考資料をVisual C++ 2005 Express Edition版に変更

著者プロフィール

  • 石立 喬(イシダテ タカシ)

    1955年東京工大卒。同年、NECへ入社し、NEC初のコンピュータの開発に参画。磁気メモリ、半導体メモリの開発、LSI設計などを経て、1989年帝京大学理工学部教授。情報、通信、電子関係の教育を担当。2002年定年により退職し現在に至る。2000年より、Webサイト「Visual C++の勉強部屋」...

バックナンバー

連載:Javaで学ぶグラフィックス処理

もっと読む

All contents copyright © 2005-2017 Shoeisha Co., Ltd. All rights reserved. ver.1.5