Tercel::Diary

ProcessingでXファイルを解析（スタティックメッシュ篇その2）

2012-02-28T18:26:00.004+09:00

昨日の続きです。ほんのちょっと進歩してマテリアルを反映できるようになりました。

こういうかっこいい戦車も……

Metasequoia にプリセットされている「Tank」

ちゃんと質感を保ったまま^※ Processing に持ってくる事ができました。やったー！

Processing で表示した結果

※ というのは大ウソで、法線無視、フラットシェーディング一択という男らしい割り切り方をしています。P3D モードにおけるレンダリングのクオリティならばこのくらいばっさり行っても許容範囲な気がします（おい

ちなみにここだけの話、X ファイル内のジオメトリの頂点数とテクスチャのUV頂点数が一致していないと表示できないという問題点があるのですが、Metasequoia が吐く X ファイルでは特に目立った問題は起きていないため、思い切って無視しています。

それでも泥縄式に拡張しただけあって、いい感じにソースが汚くなって参りました。

【Xファイル解析関連のコード】

// -------------------------------------------------
// =================================================
// Xファイル形式で記述されたメッシュを表示するテスト
//                スタティックメッシュ篇やや不完全版
// =================================================
// -------------------------------------------------


// =================================================

// メッシュ
class Mesh {
  List<Polygon> polygons;
  
  public Mesh() {
    polygons = new ArrayList<Polygon>();
  }
  
  public void add(Polygon p) {
    polygons.add(p);
  }
  
  public Polygon getPolygon(int index) {
    return polygons.get(index);
  }
  
  public void render() {
    if(polygons == null) return;
    
    for(Polygon p : polygons) p.render();
  }
}

// =================================================

// ポリゴン
class Polygon {
  List<PVector> vertices;
  List<PVector> textureCoords;
  Material      material;
  
  public Polygon() {
    vertices = new ArrayList<PVector>();
  }
  
  public void addVertex(PVector v) {
    vertices.add(v);
  }
  
  public void addTextureCoord(PVector t) {
    if(textureCoords == null) textureCoords = new ArrayList<PVector>();
    textureCoords.add(t);
  }
  
  public void setMaterial(Material m) {
    material = m;
  }
  
  public void render() {
    if(vertices == null) return;
    
    if(vertices.size() == 3)      beginShape(TRIANGLES);
    else if(vertices.size() == 4) beginShape(QUADS);
    else                          beginShape();
    
    pushStyle();
    
    // マテリアルの設定
    if(material != null) {
      noStroke();
      fill(material.faceColor);

      shininess(material.power);
      specular(material.specularColor);
      emissive(material.emissiveColor);      
      
      if(material.texture != null) {
        // テクスチャが存在する場合
        texture(material.texture);
      }
    }
    
    for(int i = 0; i < vertices.size(); ++i) {
      PVector v = vertices.get(i);
      // 実際の頂点とテクスチャのUVが一致していない場合は表示できないよ
      if(textureCoords != null && 
          textureCoords.size() == vertices.size() && 
            material != null &&
              material.texture != null) {
          
        PVector tex = textureCoords.get(i);
        textureMode(NORMALIZED);
        vertex(v.x, v.y, v.z, tex.x, tex.y);
      } else {
        vertex(v.x, v.y, v.z);
      }
    }
    popStyle();
    endShape();
  }
}

// =================================================

// マテリアル
class Material {
  int     faceColor;
  float   power;
  int     specularColor;
  int     emissiveColor;
  PImage  texture;
  
  public Material() {
    texture       = null;
  }
}

// =================================================

// Xファイルローダクラス
class Loader {
  // Xファイルの文字列データが一括して格納される
  private StringBuffer data;
  
  // Xファイル内の部分文字列の先頭を指すインデックス
  private int          index;

  // メッシュとマテリアルの連想配列
  private HashMap<String, Mesh>     meshMap;
  private HashMap<String, Material> materialMap;

  // ポリゴンごとの頂点インデックス格納用
  private List<List<Integer>> vertexIndices;

  // メッシュ
  private Mesh mesh;

  // ----------------------------------------
  public Loader() {
    data = new StringBuffer();
    
    vertexIndices = new ArrayList<List<Integer>>();
    
    meshMap      = new HashMap<String, Mesh>();
    materialMap  = new HashMap<String, Material>();
  }
  
  // ----------------------------------------
  // ファイル名を指定してメッシュを読み込み 
  public Mesh loadX(final String fileName) { 

    // 色々初期化
    mesh = null;

    vertexIndices.clear();
    meshMap.clear();
    materialMap.clear();
    
    index = 0;

    data.delete(0, data.length());
    String[] lines = loadStrings(fileName);
    for(String line : lines) {
      // コメント除去の処理
      data.append(line.replaceAll("(##|//).*", "") + "\n");
    }
    
    List<PVector> textureCoordsList = null;
    while(index < data.length()) {
      String token = getNextToken();
      
      // テンプレートは読み飛ばす
      if(token.equals("template")) skipBlock();

      // ====================
      // Material（前方宣言）
      // ====================
      else if(token.equals("Material")) loadMaterialData();

      // ====================
      // MESH
      // ====================
      else if(token.equals("Mesh")) mesh = loadMeshData();

      // ====================
      // MeshTextureCoords
      // ====================
      else if(token.equals("MeshTextureCoords")) loadTextureCoords();
      
      // ====================
      // MeshMaterialList
      // ====================
      else if(token.equals("MeshMaterialList")) loadMaterialListData();

    }
    return mesh;
  }

  // ----------------------------------------
  // { } で囲まれたブロックを読み飛ばす
  private void skipBlock() {
    try {
      // はじめの { まで読み飛ばす 
      while(!(getNextToken().equals("{")) && index < data.length()) ;
      
      int n = 1;
      while(n > 0 && index < data.length()) {
        String token = getNextToken();
        if(token.equals("{")) ++n;
        if(token.equals("}")) --n;
      }
    } catch (RuntimeException ex) {
      ex.printStackTrace();
    }
  }
  
  // ----------------------------------------
  // Meshブロックを読む
  private Mesh loadMeshData() {
    Mesh mesh = new Mesh();

    String token = getNextToken();
    // トークンが { でない場合
    if(!token.equals("{")) {
      meshMap.put(token, mesh);
      getNextToken();
    }
        
    // 頂点リストを作成
    List<PVector> verticesList = new ArrayList<PVector>();
    vertexIndices.clear();
        
    // 全頂点数を取得
    int nVertices = Integer.parseInt(getNextToken());
        
    // 各頂点の座標を取得
    for(int i = 0; i < nVertices; ++i) {
      float x = Float.parseFloat(getNextToken());
      float y = Float.parseFloat(getNextToken());
      float z = Float.parseFloat(getNextToken());          
      verticesList.add(new PVector(x, y, z));
    }
        
    // ポリゴン数を取得
    int nPolygons = Integer.parseInt(getNextToken());
    for(int i = 0; i < nPolygons; ++i) {
      Polygon poly = new Polygon();
      
      // ポリゴンごとに頂点インデックスリストを作る
      List<Integer> indices = new ArrayList<Integer>();
      
      // 1ポリゴンあたりの頂点数
      int nVerticesPerPolygon = Integer.parseInt(getNextToken());
      for(int j = 0; j < nVerticesPerPolygon; ++j) {
            
        // 頂点インデックスを取得
        int verticesIndex = Integer.parseInt(getNextToken());

        indices.add(verticesIndex);
        
        // 頂点インデックスから該当する頂点を取得してポリゴンに追加
        poly.addVertex(verticesList.get(verticesIndex));
      }
      
      vertexIndices.add(indices);
      mesh.add(poly);
    }
    return mesh;
  }


  // ----------------------------------------
  // Materialブロックを読む
  private Material loadMaterialData() {
    Material material = new Material();

    String token = getNextToken();
    // トークンが { でない場合
    if(!token.equals("{")) {
      materialMap.put(token, material);
      getNextToken();
    }
    
    int r = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    int g = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    int b = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    int a = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    
    material.faceColor = a << 24 | r << 16 | g << 8 | b;
    
    material.power = Float.parseFloat(getNextToken());

    r = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    g = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    b = 0xFF & round(255 * Float.parseFloat(getNextToken()));

    material.specularColor = 0xFF << 24 | r << 16 | g << 8 | b;

    r = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    g = 0xFF & round(255 * Float.parseFloat(getNextToken()));
    b = 0xFF & round(255 * Float.parseFloat(getNextToken()));

    material.emissiveColor = 0xFF << 24 | r << 16 | g << 8 | b;

    if(getNextToken().equals("TextureFilename")) {
      
      getNextToken();  // トークン「{」
      String textureFileName = getNextToken();
      material.texture = loadImage(textureFileName);
    }
    return material;
  }

  // ----------------------------------------
  // MeshTextureCoordsブロックを読む
  private void loadTextureCoords() {
    getNextToken();  // 「{」
    int nCoords = Integer.parseInt(getNextToken());
    List<PVector> coordsList = new ArrayList<PVector>();
    
    if(mesh == null) return;
    
    for(int i = 0; i < nCoords; ++i) {
      float u = Float.parseFloat(getNextToken());
      float v = Float.parseFloat(getNextToken());
      
      PVector texCoords = new PVector(u, v);
      coordsList.add(texCoords);
    }
    
    
    for(int i = 0; i < vertexIndices.size(); ++i) {
      Polygon       poly    = mesh.getPolygon(i);
      List<Integer> indices = vertexIndices.get(i);
      for(Integer index : indices) {
        poly.addTextureCoord(coordsList.get(index));
      }
    }
    
  }

  // ----------------------------------------
  // MeshMaterialListブロックを読む
  private void loadMaterialListData() {
    getNextToken();  // 「{」
    
    int nMaterials = Integer.parseInt(getNextToken());
    int nPolygons  = Integer.parseInt(getNextToken());
    
    // マテリアルインデックスリストを作成
    List<Integer> materialIndicesList = new ArrayList<Integer>();
    for(int i = 0; i < nPolygons; ++i)
      materialIndicesList.add(Integer.parseInt(getNextToken()));

    // Xで定義されているマテリアルを読み込み    
    List<Material> materialList = new ArrayList<Material>();
    for(int i = 0; i < nMaterials; ++i) {
      String token = getNextToken();
      while(token.equals("{")) token = getNextToken();
      while(token.equals("}")) token = getNextToken();
      
      Material material;
      if(token.equals("Material")) material = loadMaterialData();
      else                         material = materialMap.get(token);
      materialList.add(material);
    }
    
    for(int i = 0; i < nPolygons; ++i) {
      Polygon poly = mesh.getPolygon(i);
      int index    = materialIndicesList.get(i);
      Material m   = materialList.get(index);
      
      poly.setMaterial(m);
    }
    
  }

  // ----------------------------------------
  // 次のトークンを得る
  private String getNextToken() {
    // 最初にインデックスの範囲チェック
    if(!(index < data.length())) return "";  // ファイルの末尾に達した場合
    
    // デリミタ（区切り文字）の読み飛ばし
    String delimiters = " \t\r\n,;\"";
    while(index < data.length() && 
          strChr(delimiters, data.charAt(index)))  
      ++index;

    if(!(index < data.length())) return "";  // ファイルの末尾
 
    // 中括弧の検出
    char[] c = { data.charAt(index) };
    if(c[0] == '{' || c[0] == '}') {
      ++index;
      return new String(c);
    }
    
    // それ以外のトークン検出
    delimiters += "{}";
    int startIndex = index;
    int endIndex   = index;
    while(endIndex < data.length() &&
          !strChr(delimiters, data.charAt(endIndex)))
      ++endIndex;
    index = endIndex;
    return data.substring(startIndex, endIndex);
  }
  
  // ----------------------------------------
  // 文字列シーケンスsの中に、cが含まれているかどうかをチェック
  private boolean strChr(final CharSequence s, final char c) {
    for(int i = 0; i < s.length(); ++i)
      if(s.charAt(i) == c) return true;
    return false;
  }
}

【使い方】

Mesh mesh;

void setup() {
  size(640, 480, P3D);
  
  // ファイル名を指定してXをロード  
  mesh = new Loader().loadX("3 テクスチャ付き.x");
  
}

void draw() {
  background(0);
  camera();
  lights();
  noStroke();

  camera(-500, 0, -500, 0, 0, 0, 0, -1, 0);
  scale(200);
  
  // メッシュの描画
  mesh.render();
}

ちなみにエラー処理とかはあまりしていないので、X ファイルが規格外だったりテクスチャファイルが見つからなかったりすると平気で RuntimeException をぶん投げると思いますがそこはご愛嬌という事で。

【MESH GURUのサンプルに挑戦】

鎌田茂雄著『MESH GURU』の第 8 章で使用されているサンプルのうち、アニメーションを含まない X ファイルをパースしてみる事にします。以下では、Processing における実行結果と、使用した X ファイルをそのまま掲載しています。

『MESH GURU』のサンプルは Metasequoia が吐く X ファイルとはじゃっかん構造が違い、冒頭で名前つきの Material を定義して、以降は識別子で Material を参照する形になっています。なんとも意地の悪い書き方ではありますが、可能な限り対応する事にしました。

【世界最小Xファイル】

xof 0303txt 0032

///////////////////メッシュ//////////////////////////
Mesh Mesh_Triangle
{
  //////頂点データ部////////////
   3;
  -1.0;-1.0;0.0;,
  -1.0;1.0;0.0;,
  1.0;-1.0;0.0;;
  //////ポリゴンデータ部///////
  1;
  3;0,1,2;; 
}

【マテリアル付き】

xof 0303txt 0032
//////////////////////////マテリアル/////////////////////
Material Triangle_Blue
{
  0.0;0.0;1.0;1.0;;
  51.2;
  0.0;0.0;0.0;;
  0.0;0.0;0.0;;
}
///////////////////メッシュ//////////////////////////
Mesh Mesh_Triangle
{
   3;
  -1.0;-1.0;0.0;,
  -1.0;1.0;0.0;,
  1.0;-1.0;0.0;;
  1;
  3;0,1,2;;
  ///////////////////マテリアルリスト////////////
  MeshMaterialList 
  {
    1;
    1;
    0;
    { Triangle_Blue }
  }
}

【テクスチャ付き】

xof 0303txt 0032
//////////////////////////マテリアル/////////////////////
Material Triangle_Tex
{
  1.0;1.0;1.0;1.0;;
  51.2;
  0.0;0.0;0.0;;
  0.0;0.0;0.0;;
  TextureFilename
  {
    "コンクリート.bmp";
  }
}
///////////////////メッシュ//////////////////////////
Mesh Mesh_Triangle
{
   3;
  -1.0;-1.0;0.0;,
  -1.0;1.0;0.0;,
  1.0;-1.0;0.0;;
  1;
  3;0,1,2;;
  ///////////////////テクスチャ座標///////////////
  MeshTextureCoords
  {
    3;
    0.0;1.0;,
    0.0;0.0;,
    1.0;1.0;;
  }
  ///////////////////マテリアルリスト////////////
  MeshMaterialList 
  {
    1;
    1;
    0;
    { Triangle_Tex }
  }
}

※ 文字コードが UTF-8 以外の場合、日本語のテクスチャを読み込もうとした瞬間に死にます。

【おまけ：魔道少女】

↓

【おまけその2：どせいさん】

こちらのサイトから拝借した「どせいさん」を表示してみました。なお、左右が反転しているのは、座標軸の取り方によるものです。

ProcessingでXファイルを解析（スタティックメッシュ篇その1）

2012-02-27T21:04:00.000+09:00

今日は、DirectX で標準的に使用されていた X ファイルを読み込んで、Processing で表示してみたいと思います。すなわち ——

Metasequoia にプリセットされている「日本橋麒麟像」

こんなふうに任意のモデリングソフトで作成したデータを、

上図のように、Processing でそのまま使えるようにしてしまおうという試みです。

しかし僕がアホすぎるせいで、まだ静止した 3D モデルの幾何形状だけしか読めていませんごめんなさい（アニメーションはおろか物体の色付けやテクスチャも無し）。

この課題の一番のネックは X ファイルのパーサです。

X のフォーマットはそれなりに自由度が高く、使用するモデリングソフトによって構造が変わったりします。すべての状況に適応できる汎用性の高いコードを本気で書こうとすると、割と本格的な構文解析器を実装する必要があるような気がして心が折れます。

というわけで、個人的にはひとまず姑息的なパーサを手っ取り早く作り、なんか問題が起きたら適宜設計を見直しつつ修正していこうかと思います。

とりあえずXファイル読み込み関連をまとめたプログラム（長いので折りたたんであります）。

// -------------------------------------------------
// =================================================
// Xファイル形式で記述されたメッシュを表示するテスト
//                    スタティックメッシュ篇不完全版
// =================================================
// -------------------------------------------------


// =================================================

// メッシュ
class Mesh {
  List<Polygon> polygons;
  
  public Mesh() {
    polygons = new ArrayList<Polygon>();
  }
  
  public void add(Polygon p) {
    polygons.add(p);
  }
  
  public void render() {
    if(polygons == null) return;
    
    for(Polygon p : polygons) p.render();
  }
}

// =================================================

// ポリゴン
class Polygon {
  List<PVector> vertices;
  
  public Polygon() {
    vertices = new ArrayList<PVector>();
  }
  
  public void add(PVector v) {
    vertices.add(v);
  }
  
  public void render() {
    if(vertices == null) return;
    
    if(vertices.size() == 3)      beginShape(TRIANGLES);
    else if(vertices.size() == 4) beginShape(QUADS);
    else                          beginShape();
    for(PVector v : vertices) vertex(v.x, v.y, v.z);
    endShape();
  }
}

// =================================================

// Xファイルローダクラス
class Loader {
  // Xファイルの文字列データが一括して格納される
  private StringBuffer data;
  
  // Xファイル内の部分文字列の先頭を指すインデックス
  private int          index;

  // 頂点リスト
  List<PVector> verticesList;
  
  // ----------------------------------------
  public Loader() {
    data = new StringBuffer();
    verticesList = new ArrayList<PVector>();
  }
  
  // ----------------------------------------
  // ファイル名を指定してメッシュをロード 
  public Mesh loadX(final String fileName) { 
    Mesh mesh = new Mesh();
    
    index = 0;
    data.delete(0, data.length());
    String[] lines = loadStrings(fileName);
    for(String line : lines) {
      // ここで本来はコメント除去の処理が必要
      
      data.append(line + "\n");
    }
    while(index < data.length()) {
      String token = getNextToken();
      
      // テンプレートは読み飛ばす
      if(token.equals("template")) skipBlock();

      // ====================
      // MESH
      // ====================
      else if(token.equals("Mesh")) {
        token = getNextToken();
        // トークンが { でない場合
        if(!token.equals("{")) {
          getNextToken();
        } 
        /*
        else {
          // 名前つきマップに入れたい
        }*/
        
        // 頂点リストを作成
        verticesList.clear();
        
        // 全頂点数を取得
        int nVertices = Integer.parseInt(getNextToken());
        
        // 各頂点の座標を取得
        for(int i = 0; i < nVertices; ++i) {
          float x = Float.parseFloat(getNextToken());
          float y = Float.parseFloat(getNextToken());
          float z = Float.parseFloat(getNextToken());          
          verticesList.add(new PVector(x, y, z));
        }
        
        // ポリゴン数を取得
        int nPolygons = Integer.parseInt(getNextToken());
        for(int i = 0; i < nPolygons; ++i) {
          Polygon poly = new Polygon();
          
          // 1ポリゴンあたりの頂点数
          int nVerticesPerPolygon = Integer.parseInt(getNextToken());
          for(int j = 0; j < nVerticesPerPolygon; ++j) {
            
            // 頂点インデックスを取得
            int verticesIndex = Integer.parseInt(getNextToken());
            
            // 頂点インデックスから該当する頂点を取得してポリゴンに追加
            poly.add(verticesList.get(verticesIndex));
          }

          mesh.add(poly);
        }
      }
    }
    return mesh;
  }
  
  // ----------------------------------------
  // { } で囲まれたブロックを読み飛ばす
  private void skipBlock() {
    try {
      // はじめの { まで読み飛ばす 
      while(!(getNextToken().equals("{")) && index < data.length()) ;
      
      int n = 1;
      while(n > 0 && index < data.length()) {
        String token = getNextToken();
        if(token.equals("{")) ++n;
        if(token.equals("}")) --n;
      }
    } catch (RuntimeException ex) {
      ex.printStackTrace();
    }
  }
  
  // ----------------------------------------
  // 次のトークンを得る
  private String getNextToken() {
    // 最初にインデックスの範囲チェック
    if(!(index < data.length())) return "";  // ファイルの末尾に達した場合は""を返す
    
    // デリミタ（区切り文字）の読み飛ばし
    String delimiters = " \t\r\n,;\"";
    while(index < data.length() && 
          strChr(delimiters, data.charAt(index)))  
      ++index;

    if(!(index < data.length())) return "";  // ファイルの末尾
 
    // 中括弧の検出
    char[] c = { data.charAt(index) };
    if(c[0] == '{' || c[0] == '}') {
      ++index;
      return new String(c);
    }
    
    // それ以外のトークン検出
    delimiters += "{}";
    int startIndex = index;
    int endIndex   = index;
    while(endIndex < data.length() &&
          !strChr(delimiters, data.charAt(endIndex)))
      ++endIndex;
    index = endIndex;
    return data.substring(startIndex, endIndex);
  }
  
  // ----------------------------------------
  // 文字列シーケンスsの中に、cが含まれているかどうかをチェック
  private boolean strChr(final CharSequence s, final char c) {
    for(int i = 0; i < s.length(); ++i)
      if(s.charAt(i) == c) return true;
    return false;
  }
}

で、使い方はこんな感じ。

Mesh mesh;

void setup() {
  size(640, 480, P3D);
  
  // ファイル名を指定してXをロード  
  mesh = new Loader().loadX("nihonbasikirin.x");
  
}

void draw() {
  background(0);
  camera();
  lights();
  noStroke();

  camera(500, 500, 500, 0, 300, 0, 0, -1, 0);
    
  scale(300);
  
  // メッシュの描画
  mesh.render();
}

ちなみに、現在の制約では、一つの X ファイルに複数のメッシュ情報が定義されている場合、うまく読み込む事ができません（対応は簡単ですが、それが必要な例を見た事がないのでとりあえずこのままでいいや的な……）。

ただ、視覚的には複数のように見えるメッシュであっても、内部的には一つのメッシュにまとめられているような場合は、正しくロードする事ができます（下図参照）。

Metasequoia にプリセットされている「魔道少女」

これを Processing で表示させると、ちゃんと 2 つのモデルが表示されます。

これからいろいろ忙しくなりますが、暇を見つけて『色付け』と『アニメーション』に対応していけたらなぁ……と考えています。うーん……。

おもちゃであそぶ

2012-02-24T22:14:00.000+09:00

今日はおとなしくおうちで遊ぶ事にしました。

Processingでスキンメッシュアニメーション

2012-02-23T00:21:00.001+09:00

【とにかく、ぐにゃぐにゃ】

今まで黙っていましたが、実はこれまで Processing で表示させていた立体オブジェクトは、そのほとんどが剛体でした（オーロラは例外）。

というわけで、ブログをご覧の皆さんに僕が剛体しか出せない事がバレる前に、スキンメッシュという謎の技術を使って時々刻々と変形する形状の表示に挑戦してみたいと思います。

……と見せかけて、僕は根性無しなので、小規模のサンプルを作ったところで力尽きました⊂⌒~⊃｡Д｡)⊃

たいへん気色悪い物体がぐにゃぐにゃ動いているのがお解りいただけたかと思います。

スキンメッシュは、その名の通り表皮っぽいメッシュの事で、ボーンと呼ばれる仮想的な剛体の動きに追従して頂点が変形します（→OpenProcessing）。

この記事でも手法を詳しく紹介しようかと思いましたが、理論自体は大した事ないくせに解説がやたら面倒なのでやめました。理論の詳細は『ゲームエンジンアーキテクチャ』の第11章や『MESH GURU with Direct3D 10/11』の第12章あたりにざっくりまとめられていますし、Web にも良質なドキュメントがあるのでそちらを参照してください（逃

ただ、OpenProcessing で “Skinning” や “Skin Mesh” で検索してもそれっぽい先行作品がヒットしなかったので、恐らく今回のスキニングは OpenProcessing 史上初の試みだと思います。えっへん。

Processingで魔法使いの組分け帽子を作る

2012-02-20T23:44:00.000+09:00

今日はちょっとしたパターン認識をやってみようと思います。あ、タイトルは釣りですよ。

とりあえずデモ動画をご覧ください。

こんな感じで、ある「特徴」を手掛かりにして正体をつきとめる処理は「パターン認識」と呼ばれており、文字認識などに応用されている技術です。

現在に至るまで、様々な認識手法（というか識別手法）が提案されていますが、今日はフィードフォワード型のニューラルネットワーク（3層パーセプトロン）にフォーカスを当てたいと思います。

【ニューラルネットワークの構成】

ここでは、ニューラルネットワークの理論的背景について、石井健一郎ほか著『わかりやすいパターン認識』の 3 章に準拠して解説したいと思います。前提知識として、「パターン認識とは何か」という根本的な理解と、単純パーセプトロンを実装できる程度の能力があるとよいと思います。

まず、ニューロンと呼ばれる神経細胞の働きを、以下に示す閾値論理ユニットでモデル化します。これは、適当なウェイト \( \vect{w} = \left( w_0,\,\cdots,\,w_d\right)^t\) によって重みづけした入力信号 \(\vect{x} = \left(x_0,\,\cdots,\,x_d\right)^t\) の総和を基に、閾値処理を行う単純な仕組みです。出力 \(g_i(\vect{x})\) は 0 か 1 となります。

閾値論理ユニット

ニューラルネットワークは、複数の閾値論理ユニット（以下、ユニット）を含む層を、入力層から出力層まで多数並べたネットワークの事です。各ユニットは隣接する層とだけで結合しており、かつ入力層から出力層へ向かう一方向のネットワーク構成です（下図参照）。

ニューラルネットワークのある層における \(j\) 番目のユニットに注目し、これをユニット \(j\) と呼ぶ事にします。また、ユニット \(j\) の1階層前の \(i\) 番目のユニットをユニット \(i\)、1階層後ろの \(k\) 番目のユニットをユニット \(k\) と呼び、さらにユニット \(i\) からユニット \(j\) への結合の重みを \(w_{ij}\)、ユニット \(j\) からユニット \(k\) への結合の重みを \(w_{jk} \) で表します。

ここで、ネットワークにある学習パターン \(\vect{x}_p\) を入力したとき、ユニット \(i\) からの出力を \(g_{ip}\)、ユニット \(j\) への入力を \(h_{jp}\) とすると、\(h_{jp}\) は \(j\) に結合している1階層前の層内のユニットからの出力の重み付き信号の総和、つまり
\begin{align}
h_{jp}=\sum_{i}w_{ij}g_{ip}
\end{align}となります。

ユニット \(j\) の出力 \(g_{jp}\) は閾値関数 \(f_j\) を用いて
\begin{align}
g_{jp}=f_j(h_{jp})
\end{align}と表すものとします。

【誤差の定式化と重みの更新】

ニューラルネットワークを構築したら、理想的な識別結果が既知の入力データを用いて学習（重みの調整）を行う必要があります。

ある学習パターン \(\vect{x}_p\) に対する二乗誤差は、出力層のユニット \(l\) の出力 \( g_{lp} \) と、ユニット \(l\) に対する教師信号（理想的な識別結果）を \(b_{lp}\) を用いて、以下のように表す事ができます。
\begin{align}
J_p = \frac{1}{2}\sum_l \left( g_{lp} - b_{lp} \right)^2
\end{align}従って、全学習パターンに対する二乗誤差 \(J\) は
\begin{align}
J=\sum_p J_p
\end{align}となります。

識別性能を高めるためには、この誤差 \(J\) を最小化するようにユニット間の重みを調整する必要があります。つまり、ある重み \(w_{ij}\) に対して、以下の更新式を適用する事になります。
\begin{align}
w'_{ij} &= w_{ij} - \rho \frac{\partial J_p}{\partial w_{ij}} \\
&= w_{ij} - \rho \underbrace{ \frac{\partial J_p}{\partial h_{jp}}}_{\varepsilon_{jp}} \cdot \underbrace{\frac{\partial h_{jp}}{\partial w_{ij}}}_{g_{ip}} \\
&= w_{ij} - \rho \varepsilon_{jp}g_{ip}
\end{align}ここで \(\rho\) を学習係数と呼び、任意の正定数とします。

また、上式右辺の \(\varepsilon_{jp}\) は、合成関数の微分を用いて以下のように式変形できます。
\begin{align}
\varepsilon_{jp} &= \frac{\partial J_p}{\partial h_{jp}} = \frac{\partial J_p}{\partial g_{jp}} \cdot \frac{\partial g_{jp}}{\partial h_{jp}} \\
&= \frac{\partial J_p}{\partial g_{jp}} \cdot f'_j(h_{jp})
\end{align}なお、上式中の \( \dfrac{\partial J_p}{\partial g_{jp}} \) は、ユニット \(j\) が中間層か出力層かで以下のように場合分けされます。
\begin{align}
\frac{\partial J_p}{\partial g_{jp}} =
\begin{cases}
g_{jp} - b_{jp} & \textrm{（出力層）} \\
\sum_{k}\varepsilon_{kp}w_{jk} & \textrm{（中間層）}
\end{cases}
\end{align}

【閾値関数の選定】

最後に、閾値関数 \(f_j\) を考えます。ここまでの議論で、この関数が微分可能でなければならない事にはお気付きでしょう。

通常、\(f_j\) には以下のようなシグモイド関数 \(S(u)\) が選ばれます。
\begin{align}
S(u) = \frac{1}{1 + \exp(-u)}
\end{align}
この関数は、微分すると \(S'(u) = S(u)\left( 1 - S(u)\right)\) になりますから、\(f'_j\) は
\begin{align}
f'_j (h_{jp}) = g_{jp} (1 - g_{jp})
\end{align}となります。

以上を整理すると、\(\varepsilon_{jp}\) に関する再帰式を以下のように得る事ができます。
\begin{align}
\varepsilon_{jp} =
\begin{cases}
\left( g_{jp} - b_{jp} \right) g_{jp} \left( 1 - g_{jp} \right) & \textrm{（出力層）} \\
\left( \sum_k \varepsilon_{kp} w_{jk} \right) g_{jp} \left( 1 - g_{jp} \right) & \textrm{（中間層）}
\end{cases}
\end{align}上式において、\(0 < g_{jp} \left( 1 - g_{jp} \right) < 1\) であり、ユニットの出力値 \(g_{jp} \) が 0.5 のときに重みの修正量が最も大きく、0 か 1 に近づくほど修正量が小さくなります。

学習パターンが入力されると出力層で各ユニットの出力と教師信号の誤差が計算され、それに基づいて \(\varepsilon_{jp}\) が求められ、出力層での重みの修正が行われます。

この結果によって、1つ前の階層でも重みの修正が行われます。これを同様に繰り返す事で、すべての層で重みが修正されます。このように、出力層から入力層に向かって誤差を伝播させる学習こそが、誤差逆伝播法と呼ばれる所以です。

【実装】

3層から成るフィードフォワード型ニューラルネットワークの実装は以下の通りです。こちらも、David M.Bourg, Glenn Seemann 著『ゲーム開発者のための AI 入門』を参考にしました（そこそこ長いので展開してご覧ください）。

class NeuralNetworkLayer {
  int        numberOfNodes;       // 当該レイヤのノード数
  int        numberOfChildNodes;  // 子レイヤのノード数
  int        numberOfParentNodes; // 親レイヤのノード数
  
  double[][] weights;             // 親子レイヤを接続するノードの重み値
  double[][] weightChanges;       // 重み値の調整に使われる値
  
  double[]   neuronValues;        // 計算によって得られたニューロンの値
  double[]   desiredValues;       // 目的となるニューロンの値
  double[]   errors;              // 各ニューロンに関連するエラー値
  double[]   biasWeights;         // 各ニューロンに接続するバイアス重み
  double[]   biasValues;          // バイアス値 通常は +1 または -1
  double     learningRate;        // 重みの調整を計算する学習率
  
  boolean    linearOutput;        // 線型活性化関数を使用するかどうか
                                  // false ならロジスティック関数
  boolean    useMomentum;         // 重みの調整にモーメンタムを使用するかどうか
                                  // デフォルトはfalse
  double     momentumFactor;      // モーメンタム因数
  
  NeuralNetworkLayer parentLayer; // 親レイヤへの参照 入力層の場合は null
  NeuralNetworkLayer childLayer;  // 子レイヤへの参照 出力層の場合は null
  
  // コンストラクタ
  NeuralNetworkLayer() {
    parentLayer    = null;
    childLayer     = null;
    linearOutput   = false;
    momentumFactor = 0.9;
  }
  
  // 初期化メソッド
  void initialize( // int numNodes, 
                  NeuralNetworkLayer parent, 
                  NeuralNetworkLayer child) {

    neuronValues  = new double[numberOfNodes];
    desiredValues = new double[numberOfNodes];
    errors        = new double[numberOfNodes];
    
    if(parent != null) {
      parentLayer = parent;
    }
    
    if(child != null) {
      childLayer = child;
      
      weights       = new double[numberOfNodes][numberOfChildNodes];
      weightChanges = new double[numberOfNodes][numberOfChildNodes];

      biasValues = new double[numberOfChildNodes];
      biasWeights = new double[numberOfChildNodes];
      
    }
    
    // すべてに0が含まれるようにする
    for(int i = 0; i < numberOfNodes; i++) {
      neuronValues[i] = 0;
      desiredValues[i] = 0;
      errors[i] = 0;
      
      if(childLayer != null)
        for(int j = 0; j < numberOfChildNodes; j++) {
          weights[i][j] = 0;
          weightChanges[i][j] = 0;
        }
    }
    
    // バイアス値と重みを初期化する
    if(childLayer != null)
      for(int j = 0; j < numberOfChildNodes; j++) {
        biasValues[j]  = -1;
        biasWeights[j] = 0;
      }
  }
  
  void randomizeWeights() {
    for(int i = 0; i < numberOfNodes; i++) {
      for(int j = 0; j < numberOfChildNodes; j++) {
        weights[i][j] = random(-1.0, 1.0);
      }
    }
    
    for(int j = 0; j < numberOfChildNodes; j++) {
      biasWeights[j] = random(-1.0, 1.0);
    }   
  }
  
  void calculateNeuronValues() {
    if(parentLayer != null) {
      for(int j = 0; j < numberOfNodes; j++) {
        double x = 0;
        for(int i = 0; i < numberOfParentNodes; i++) {
          x += parentLayer.neuronValues[i] * parentLayer.weights[i][j];
        }
        x += parentLayer.biasValues[j] * parentLayer.biasWeights[j];
        
        if((childLayer == null) && linearOutput)
          neuronValues[j] = x;
        else
          neuronValues[j] = 1.0 / (1.0 + Math.exp(-x));
      }
    }
  }
  
  void calculateErrors() {
    if(childLayer == null) {
      // 出力層
      for(int i = 0; i < numberOfNodes; i++) {
        errors[i] = (desiredValues[i] - neuronValues[i]) *
                    neuronValues[i] * (1.0f - neuronValues[i]);
      }
    } else if (parentLayer == null) {
      // 入力層
      for(int i = 0; i < numberOfNodes; i++) {
        errors[i] = 0.0f;
      }
    } else {
      // 隠れ層
      for(int i = 0; i < numberOfNodes; i++) {
        double sum = 0;
        for(int j = 0; j < numberOfChildNodes; j++) {
          sum += childLayer.errors[j] * weights[i][j];
        }
        errors[i] = sum * neuronValues[i] * (1.0f - neuronValues[i]);
      }
    }
  }
  
  void adjustWeights() {
    if(childLayer != null) {
      for(int i = 0; i < numberOfNodes; i++) {
        for(int j = 0; j < numberOfChildNodes; j++) {
          double dw = learningRate * childLayer.errors[j] * neuronValues[i];
          if(useMomentum) {
            weights[i][j] += dw + momentumFactor * weightChanges[i][j];
            weightChanges[i][j] = dw;
          } else {
            weights[i][j] += dw;
          }
        }
      }
      
      for(int j = 0; j < numberOfChildNodes; j++) {
        biasWeights[j] += learningRate * childLayer.errors[j] * biasValues[j];
      }
    }
  }
}

class NeuralNetwork {
  NeuralNetworkLayer inputLayer;
  NeuralNetworkLayer hiddenLayer;
  NeuralNetworkLayer outputLayer;
  
  void initialize(int nNodesInput, int nNodesHidden, int nNodesOutput) {
    inputLayer  = new NeuralNetworkLayer();
    hiddenLayer = new NeuralNetworkLayer();
    outputLayer = new NeuralNetworkLayer();
    
    inputLayer.numberOfNodes       = nNodesInput;
    inputLayer.numberOfChildNodes  = nNodesHidden;
    inputLayer.numberOfParentNodes = 0;
    inputLayer.initialize(null, hiddenLayer);
    inputLayer.randomizeWeights();
    
    hiddenLayer.numberOfNodes       = nNodesHidden;
    hiddenLayer.numberOfChildNodes  = nNodesOutput;
    hiddenLayer.numberOfParentNodes = nNodesInput;
    hiddenLayer.initialize(inputLayer, outputLayer);
    hiddenLayer.randomizeWeights();
    
    outputLayer.numberOfNodes       = nNodesOutput;
    outputLayer.numberOfChildNodes  = 0;
    outputLayer.numberOfParentNodes = nNodesHidden;
    outputLayer.initialize(hiddenLayer, null);
  }
  
  void setInput(int i, double value) {
    if((i >= 0) && ( i < inputLayer.numberOfNodes)) {
      inputLayer.neuronValues[i] = value;
    }
  }
  
  double getOutput(int i) {
    if((i >= 0) && (i < outputLayer.numberOfNodes)) {
      return outputLayer.neuronValues[i];
    }
    
    return Double.MAX_VALUE;  // エラーを示す
  }
  
  void setDesiredOutput(int i, double value) {
    if((i >= 0) && (i < outputLayer.numberOfNodes)) {
      outputLayer.desiredValues[i] = value;
    }
  }
  
  void feedForward() {
    inputLayer.calculateNeuronValues();
    hiddenLayer.calculateNeuronValues();
    outputLayer.calculateNeuronValues();
  }
  
  void backPropagate() {
    outputLayer.calculateErrors();
    hiddenLayer.calculateErrors();
    
    hiddenLayer.adjustWeights();
    inputLayer.adjustWeights();
  }
  
  int getMaxOutputID() {
    double maxval = outputLayer.neuronValues[0];
    int id = 0;
    for(int i = 1; i < outputLayer.numberOfNodes; i++) {
      if(outputLayer.neuronValues[i] > maxval) {
        maxval = outputLayer.neuronValues[i];
        id = i;
      }
    }
    return id;
  }
  
  double calculateError() {
    double error = 0;
    
    for(int i = 0; i < outputLayer.numberOfNodes; i++) {
      error += Math.pow(outputLayer.neuronValues[i] - outputLayer.desiredValues[i], 2);
    }
    
    error = error / outputLayer.numberOfNodes;
    return error;
  }
  
  void setLearningRate(double rate) {
    inputLayer.learningRate  = rate;
    hiddenLayer.learningRate = rate;
    outputLayer.learningRate = rate;
  }
  
  void setLinearOutput(boolean useLinear) {
    inputLayer.linearOutput  = useLinear;
    hiddenLayer.linearOutput = useLinear;
    outputLayer.linearOutput = useLinear;
  }
  
  void setMomentum(boolean useMomentum, double factor) {
    inputLayer.useMomentum  = useMomentum;
    hiddenLayer.useMomentum = useMomentum;
    outputLayer.useMomentum = useMomentum;
    
    inputLayer.momentumFactor  = factor;
    hiddenLayer.momentumFactor = factor;
    outputLayer.momentumFactor = factor;
  }
  
  void dumpData(String filename) {
    List<String> data = new ArrayList<String>();
    data.add("--------------------------------------------------------");
    data.add("Input Layer");
    data.add("--------------------------------------------------------");
    data.add("");
    data.add("Node Values:");
    for(int i = 0; i < inputLayer.numberOfNodes; i++)
      data.add("(" + i + ") = " + inputLayer.neuronValues[i]);
    data.add("");
    data.add("Weights:");
    data.add("");
    for(int i = 0; i < inputLayer.numberOfNodes; i++)
      for(int j = 0; j < inputLayer.numberOfChildNodes; j++) 
        data.add("(" + i + ", " + j + ") = " + inputLayer.weights[i][j]);
    data.add("");
    data.add("Bias Weights:");
    for(int j = 0; j < inputLayer.numberOfChildNodes; j++)
      data.add("(" + j + ") = " + inputLayer.biasWeights[j]);
    
    data.add("");
    data.add("");
    
    data.add("--------------------------------------------------------");
    data.add("Hidden Layer");
    data.add("--------------------------------------------------------");
    data.add("");
    data.add("Weights:");
    data.add("");
    for(int i = 0; i < hiddenLayer.numberOfNodes; i++)
      for(int j = 0; j < hiddenLayer.numberOfChildNodes; j++) 
        data.add("(" + i + ", " + j + ") = " + hiddenLayer.weights[i][j]);
    data.add("");
    data.add("Bias Weights:");
    for(int j = 0; j < hiddenLayer.numberOfChildNodes; j++)
      data.add("(" + j + ") = " + hiddenLayer.biasWeights[j]);

    data.add("");
    data.add("");
    
    data.add("--------------------------------------------------------");
    data.add("Output Layer");
    data.add("--------------------------------------------------------");
    data.add("");
    data.add("Node Values:");
    for(int i = 0; i < outputLayer.numberOfNodes; i++)
      data.add("(" + i + ") = " + outputLayer.neuronValues[i]);
    data.add("");
    
    String[] strings = new String[data.size()];
    data.toArray(strings);
    saveStrings(filename, strings);
  }
}

で、使い方はこんな感じです。

// --------------------------------------
// ======================================
// ニューラルネットワークを用いた
//                  2クラス識別のサンプル
// ======================================
// --------------------------------------

// NeuralNetworkオブジェクトを宣言
NeuralNetwork theBrain;

// 学習用データセット
List<PVector> class1;
List<PVector> class2;

void setup() {
  size(400, 300);
  
  // ======================================
  // 学習用データ（教師信号）のセットアップ
  // ======================================
  class1 = new ArrayList<PVector>();
  class2 = new ArrayList<PVector>();
  
  for(int i = 0; i < 10; i++) {
    class1.add(new PVector(random(width / 2),        random(height)));
    class2.add(new PVector(random(width / 2, width), random(height)));
  }
  
  
  // ======================================
  // ニューラルネットワークのセットアップ
  // ======================================
  theBrain = new NeuralNetwork();
  
  // 入力層、隠れ層、出力層の数をそれぞれ指定
  theBrain.initialize(2, 10, 3);
  
  // 学習係数を0.2に設定
  theBrain.setLearningRate(0.2);
  
  // 局所解回避のためのモーメンタムを使用
  theBrain.setMomentum(true, 0.9);
  
  // ======================================
  // 学習
  // ======================================
  double error = 1;
  long c = 0;

  // 収束するか、一定の回数に達するまで反復
  while((error > 0.01) && (c < 50000)) {
    error = 0;
    c++;

    // ------------------------------------
    // クラス1の学習
    // ------------------------------------
    for(int i = 0; i < class1.size(); i++) {
      // 入力層に、インデックスとデータの組を与える
      // ※ データは必ず 0 ～ 1 の範囲
      theBrain.setInput(0, (float)class1.get(i).x / width);
      theBrain.setInput(1, (float)class1.get(i).y / height);
      
      // 出力層に、インデックスとデータの組を与える
      // ※ シグモイド関数を使用しているため、
      //      0.1 → 不活性
      //      0.9 → 活性
      //    となるようにデータを与える
      theBrain.setDesiredOutput(0, 0.9);
      theBrain.setDesiredOutput(1, 0.1);
      
      // 誤差の計算
      theBrain.feedForward();
      error += theBrain.calculateError();
      
      // 誤差逆伝播
      theBrain.backPropagate();
    }
    
    // ------------------------------------
    // クラス2の学習
    // ------------------------------------
    for(int i = 0; i < class2.size(); i++) {
      theBrain.setInput(0, (float)class2.get(i).x / width);
      theBrain.setInput(1, (float)class2.get(i).y / height);
      
      theBrain.setDesiredOutput(0, 0.1);
      theBrain.setDesiredOutput(1, 0.9);
      
      theBrain.feedForward();
      error += theBrain.calculateError();
      theBrain.backPropagate();
    }
    
    // 誤差の平均化
    error = error / (float)(class1.size() + class2.size());
  }  
}


void draw() {
  background(0);
  
  // ======================================
  // 識別境界の視覚化
  // ======================================
  for(int i = 0; i < width; i++) {
    for(int j = 0; j < height; j++) {
      
      // ------------------------------------
      // データ入力と識別
      // ------------------------------------
      theBrain.setInput(0, (float)i / width);
      theBrain.setInput(1, (float)j / height);
      theBrain.feedForward();
      
      // 識別結果の取得
      int ret = theBrain.getMaxOutputID();
      
      if(ret == 0)      stroke(0xFFFF7F7F);  // クラス1の場合
      else if(ret == 1) stroke(0xFF7F7FFF);  // クラス2の場合
      
      // 識別境界の表示
      point(i, j);
    }
  }
  
  // 学習データの表示
  pushStyle();
  strokeWeight(3);
  
  stroke(0xFFFF0000);
  for(PVector c1 : class1) point(c1.x, c1.y);

  stroke(0xFF0000FF);
  for(PVector c2 : class2) point(c2.x, c2.y);

  popStyle();
}

ProcessingでAI

2012-02-18T22:13:00.002+09:00

【Boid 理論を眺める】

くろねこさんが Boid 理論を使ってパーティクルを動かしていたのが気になって、僕もちょっと背伸びをして大風呂敷を広げてみました。……あ、タイトルは釣りです。

Boid 理論とは、団体行動を取っているような「群れ」を再現するシミュレーションです。ご存知ない方は、先入観なしに以下の動画をご覧ください。

この理論の元ネタは、Craing Reynolds によって 1987 年に発表された “Flocks, Herds, and Schools: A Distributed Behavioral Model” という論文で、以下の行動パターンを組み合わせるだけで集団行動を巧妙にシミュレートできるというものです。

Cohesion （群れにくっつくよ！）
Alignment （みんなで同じ方向を向くよ！）
Separation （ぶつかりそうになったらちょっと離れるよ）

おもしろいですね。

なお、サンプル動画の中では、個体の「視野角の広さ」を切り替えています。

各個体が広い視野角を持つ場合、群れは横に大きく膨らむ形になります。

一方、視野角が狭い場合、個体は一列に並びます。

なお、今回の実装にあたって参考にした文献は、David M. Bourg, Glenn Seemann 著『ゲーム開発者のための AI 入門』です。

ほとんど移植に近いので恐縮ですが、お勉強日記なのでたまにはこういうのもいいかなーと。ちなみに実際に動作するサンプルとソースコードは OpenProcessing の方に上げてありますが、ほとんど C++ 実装からの移植なので微妙かも知れません。

ProcessingとArduinoをくっつける

2012-02-16T20:59:00.001+09:00

【本日のお便りコーナー】

twitter4Jかー。この前購入してたarduino使うのかと思ってた∩（・ω・）∩2012年2月16日 9:7 via HootSuiteくろねこ
laprasDrum

—— というわけで、今日は Processing と Arduino を合体させてこんなの（↓）を作ってみようと思います。わー（歓声）ぱちぱち（拍手） ……あ、ちなみにこれの正体は最後に明らかになりますよ。

—— とは言ったものの、Processing と Arduino を連繋させる方法が見当もつきません。どこかにいいサンプルはないかと探していたら、こちらの記事を見つけたので参考にする事にしました。

【配線】

まずは Arduino の配線です。13 番ピンに LED を挿すらしいです（※ ピンぼけ容赦）。

LEDのアノードを13番ピンに、カソードをGNDに接続

【スケッチ（クリックで拡大）】

左側が Arduino スケッチ、右側が Processing スケッチです。サンプルとほとんど一緒。

今回は、配線 → Arduinoスケッチのアップロード → Processing スケッチの実行の順に行いました。あまり厳密にやらなくてもよさそうですが。

注意点としては、Processing 側の以下の箇所。

myPort = new Serial(this, Serial.list()[1], 9600);

コンストラクタの第2引数には Arduino が接続されているポートを指定してやる必要があります。

ここで使用しているSerial クラスの list(); メソッドは、使用可能なポート名の配列を返すメソッドです。その中から適切なポートを選んでやればよいみたいですね（僕の環境では [1] ）。

もしも何番のポートを選べばよいか判らない場合は、一度 Serial.list(); を実行して有効なポート名を列挙し、Arduino が接続されているポートと一致するインデックスを探すのがよいと思います。

【実行結果】

マウスボタンを押下すると LED が点灯し、離すと消灯します。

次に、Arduino からデータを送ってみる事にします。

……で、そのサンプルがこちらなのですが、残念ながら僕は可変抵抗器を持っていません。しかも、近所に電子部品を売っているお店がないので、とりあえず CdS セルで代用してみようと思います。

まずは CdS セルのおさらいも兼ねて、Arduino だけをいじってみる事にします。

【配線】

CdS セルと LED を使ってこんな感じの配線をします。

実写だとこんな感じ。

画像では見づらいですが、CdS セルをアナログ入力用の 0 番ピンに、LED を出力の 9 番ピンにそれぞれ繋いでいます。

【Arduino スケッチ】

次に、以下のようなコードを書いて、Arduino に送ります。

#define LED 9
int val = 0;

void setup() {
  pinMode(LED, OUTPUT);
}

void loop() {
  val = analogRead(0);
  
  analogWrite(LED, val/4);
  delay(10);
}

上記のソースの 9 行目にある

val = analogRead(0);

は、アナログ入力ピンの電圧を調べる関数です。Arduino は AD コンバータによって、0 ～ 5 [V] の入力電圧を 0 ～ 1023 の数値に変換しているため、変数 val は 0 ～ 1023 の範囲をとる事になります。

【実行結果】

これで、明るさに応じて LED の光の強さが変わるものができました。

このサンプルで得た電圧のデータを Processing に送ってやれば、今までにできなかった事ができるようになるのではないでしょうか。

やってみましょう。

【配線】

とりあえず、LED はもう要らない子なのでブレッドボードから外し、以下のような配線に戻します。

【Arduino スケッチ】

int val = 0;

void setup() {
  Serial.begin(9600);
}

void loop() {
  val = analogRead(0);
  
  // valを0～255の範囲に正規化してP5に送る
  Serial.write(map(val, 0, 1023, 0, 255));  
  delay(100);
}

【Processing スケッチ】

import processing.serial.*;
Serial myPort;

int val;
 
void setup()
{
  size(400, 300);
  
  String port = Serial.list()[1];
  myPort = new Serial(this, port, 9600);
}
 
void draw()
{
  // 電圧の強さに応じて背景色が変わるよ
  background(val);
}
 
void serialEvent(Serial p){
  // Arduinoからデータを読み取り
  val = myPort.read();
}

【実行結果】

CdS セルに充分な光が当たっているときは下図左、センサを手で完全に覆い隠すと下図右のようになりました。

基本的に、一度に Arduino とやり取りできるデータ量が 8 bits なので、簡単のために val の範囲を 0 ～ 255 に正規化しています（0 ～ 1023 のデータも、一応がんばれば送る事はできますが）。

【おまけ】

最後に、電圧の変化を無駄にビジュアル化してみます。

Arduino のスケッチはそのままで、Processing の方を以下のように書き換えるだけです。

【Processing スケッチ】

import processing.serial.*;
Serial myPort;

PFont font;

int   index;
int[] data;

int   val;
long  time;

void setup()
{
  size(640, 480, P3D);
  
  font = createFont("Meiryo", 20);
  
  index = 0;
  data  = new int[width];
  for(int i = 0; i < data.length; ++i) data[i] = -1;
  
  String port = Serial.list()[1];
  myPort = new Serial(this, port, 9600);
}

void draw()
{
  background(0);
  camera(0, -0xFF,   width, 
         0, -0xFF/2, 0, 
         0,  1,      0);
         

  float angle = radians(time) * 0.1f;
  rotateY(angle);
  
  pushMatrix();
  translate(-width/2, -0xFF, 0);
  
  int prevY = 0xFF;
  for(int x = 1; x < width; ++x) {
    int v = data[(x + index) % data.length];  // 電圧: 0 ～ 255
    int y = 0xFF - v;  // y座標値

    stroke(255, v);
    if(prevY < 0xFF) line(x - 1, prevY, x, y);
    prevY = y;
  }
  popMatrix();
  
  // じめん
  stroke(0, 255, 0, 100);
  for(int i = -width/2; i <= width / 2; i += 40) {
    line(-width/2, 0, i, width/2, 0, i);
    line(i, 0, -width/2, i, 0, width/2);
  } 
  
  pushMatrix();
  camera();
  hint(DISABLE_DEPTH_TEST);
  textMode(SCREEN);
  // 1単位あたりの分解能: 5000[mV] / 255 ≒ 19.6[mV]
  text("Voltage: " + (val * 19.6) + " [mV]", 10, 10);
  
  hint(ENABLE_DEPTH_TEST);
  popMatrix();
  
  ++time;
}
 
void serialEvent(Serial p){
  // Arduinoからデータを読み取り
  val = myPort.read();
  
  index = (index + 1) % data.length;
  data[index] = val;
}

ProcessingとTwitterをくっつける

2012-02-15T23:52:00.000+09:00

くろねこさんが、昨日のサンプルをさっそく fork して下さいました（→くろねこさんの記事）。

パーティクルの動きの改良に加え、ソースコードもリファクタリングされて品質が上がっています。さすが！

さてさて、今日は Processing を Twitter に組み合わせてへんなのを作りますよー∩( ・ω・)∩

なんか、テキスト領域に Twitter のユーザ名を入れると、フォローしているおともだちが3次元空間を漂うというなんとも愉快なスケッチです。

実は、ほとんど昨日のソースコードを流用しただけの手抜き更新というのは内緒だよ。

【Processing と GUI】

本日一発目のテーマは、Processing スケッチへの GUI コンポーネントの組み込みです。とりあえず完成すると以下のようになります（少し見づらいですが、動画の下の方にある藍色の矩形がテキストフィールドです）。

Processing 公式サイトにはいくつかのサードパーティ製ライブラリが紹介されています。僕はその中から、ControlP5 というライブラリを選びました（現時点の最新バージョンは ControlP5 0.5.4）。

OpenGL レンダラの環境でテキストフィールドを使用するためのミニマルなサンプルコードは以下の通り。なお、外部ライブラリの設置方法に関しては、うえちょこさんの記事に詳しく書かれています。

【ソースコード（展開してご覧ください）】

import processing.opengl.*;
import controlP5.*;

ControlP5 controlP5;
Textfield tfUserName;  // テキストフィールド


void setup() {
  size(400, 300, OPENGL);
  controlP5 = new ControlP5(this);
  
  tfUserName = controlP5.addTextfield("userName", 10, 10, 200, 20);
  tfUserName.setFocus(true);
}


void draw() {
  hint(ENABLE_DEPTH_TEST);
  background(0);

  /* ============================== */
  /* ここに描画用のコードを書きます */
  /* ============================== */
  
  // GUIコンポーネント描画のためにカメラ位置をリセット
  camera();
  hint(DISABLE_DEPTH_TEST);
}


// テキストが入力された際にコールバックされるメソッド
public void userName(String txt) {
  if(txt.trim().length() < 1) return;

  println("「" + txt + "」が入力されました");
}

デモ動画のソースコードは……適当に作ったため、あまり見せられる代物ではありませんが、たぶん昔の僕だったら「汚くてもいいからソースコードを知りたい」と思ったでしょうから、こっそり晒しておきますね。

【おまけソースコード（展開してご覧ください）】

import processing.opengl.*;
import javax.media.opengl.*;  
import controlP5.*;

final int FIELD_SIZE = 1000;
final int FIELD_STEP =  100;

ControlP5 controlP5;
Textfield tfUserName;       // テキストフィールド

PFont  font;  
long   time;  

List<Message> messageList;  // メッセージリスト

void setup() {
  size(640, 480, OPENGL);
  
  messageList = new LinkedList<Message>();
  
  controlP5 = new ControlP5(this);
  
  tfUserName = controlP5.addTextfield("textfield", 10, 10, 200, 20);
  tfUserName.setFocus(true);
  
  font = createFont("Meiryo", 48);
}


void draw() {
  background(0);

  PGraphicsOpenGL pgl = (PGraphicsOpenGL)g;  

  // カメラの設定
  float angle = 0.5f * radians(time);      
  camera(300*cos(angle), -300*sin(angle), 300*sin(angle),   
         0,               0,              0,   
         0,               1,              0);
  
  // 地面の描画
  pushStyle();
  stroke(255, 0, 0);  
  for(int i = -FIELD_SIZE; i <= FIELD_SIZE; i += FIELD_STEP) { 
    line(          i, 0, -FIELD_SIZE,          i, 0, FIELD_SIZE);   
    line(-FIELD_SIZE, 0,           i, FIELD_SIZE, 0,          i);
  }
  popStyle();
  
  // α合成を有効にする
  GL gl = pgl.beginGL();
  gl.glEnable(GL.GL_BLEND);
  gl.glBlendFunc(GL.GL_SRC_ALPHA, GL.GL_ONE_MINUS_SRC_ALPHA);  
  gl.glDisable(GL.GL_DEPTH_TEST);
  pgl.endGL();
  
  // メッセージの描画
  for(Iterator iter = messageList.iterator(); iter.hasNext();) {
    Message msg = (Message)iter.next();
    // 更新して死んでたらリストから削除
    if(!msg.update()) iter.remove(); 
  }
  time++;
  
  // GUI コンポーネント描画のための設定
  camera();
  gl = pgl.beginGL();
  gl.glDisable(GL.GL_DEPTH_TEST);  // 深度テストを無効化
  gl.glEnable(GL.GL_BLEND);
  gl.glBlendFunc(GL.GL_SRC_ALPHA, GL.GL_ONE_MINUS_SRC_ALPHA);  
  pgl.endGL();
}


// テキストが入力された際にコールバックされるメソッド
public synchronized void textfield(String txt) {
  // 適当な初期位置を指定して、リストに新たなメッセージを追加する
  messageList.add(new Message(txt,
      new PVector(random(-100, 100), 0, random(-100, 100))));
}


// メッセージクラス
class Message {
  private final int   BASECOLOR = 0xFFFFFFFF;  // 基本色
  private final int   TTL       = 200;         // 生存時間
  private final float dY        = -0.5;        // y座標の増分

  private String  message;       // メッセージ
  private PVector pos;           // 位置
  private int     time;          // 経過時間  
  
  // コンストラクタ
  // -----------------------------
  // メッセージと初期位置を指定する
  Message(String msg, PVector pos) {
    this.message = msg;
    this.pos     = pos;
  }
  
  // update()
  // -----------------------------
  // 状態を更新し、生死状態を返す
  // true:  有効（生きている）
  // false: 無効（死んでいる）
  boolean update() {
    if(++time > TTL) return false;
    pos.y += dY;
    
    // ビルボーディング 
    pushMatrix();
    translate(pos.x, pos.y, pos.z);  // 物体を並進（ローカル座標系）  
      
    // モデル-ビュー行列を取得  
    PMatrix3D builboardMat = (PMatrix3D)g.getMatrix();  
    // 回転成分を単位行列に  
    builboardMat.m00 = builboardMat.m11 = builboardMat.m22 = 1;  
    builboardMat.m01 = builboardMat.m02 =   
    builboardMat.m10 = builboardMat.m12 =   
    builboardMat.m20 = builboardMat.m21 = 0;  
    
    resetMatrix();  
    applyMatrix(builboardMat);  

    // アルファ値を計算
    pushStyle();
    int a = 255 - 255 * time / TTL;
    fill(a << 24 | BASECOLOR & 0xFFFFFF);
    
    // フォントを設定
    textFont(font, 24);  
    textAlign(CENTER);
    
    // テキストを描画
    text(message, 0, 0);
    popStyle();    
       
    popMatrix();
    return true;
  }
}

【Processing と Twitter】

GUI は一旦おいといて(つ´∀｀)つ

今度は Twitter との連繋を考えてみます。ここでのゴールは、ひとまず Twitter でフォローしている方々のアイコンを 3次元空間中に散りばめる事です。

Twitter API を直叩きするのも面倒なので、Twitter4J という便利なライブラリをちょこっと使います。ライブラリは公式サイトあるいは GitHub から入手できます。今回は最新安定版である twitter4j-2.2.5.zip を選択しました。

今日はもう疲れ始めたので、OAuth 認証の要らないフォロー/フォロワー一覧を取得するところまでやろうかと思います。

JavaDoc などを参照しつつ、フレンド（自分がフォローしているアカウント）を調べるプログラムを Processing で書いてみました。

【ソースコード（展開してご覧ください）】

import twitter4j.*;

final String USER_NAME = "tercel_s";

Twitter twitter;

void setup() {
  twitter = new TwitterFactory().getInstance();
  
  try {
    
    // 指定したユーザがフォローしているアカウント一覧を取得
    IDs friends = twitter.getFriendsIDs(USER_NAME, -1);
    long[] friendsIDs = friends.getIDs();
    
    // Twitter ID からユーザ情報を逐次取得して表示
    for(long id : friendsIDs) {
      User user = twitter.showUser(id);
      println(user.getName() + " (" + user.getScreenName() + ")");
    }
    
  } catch (TwitterException ex) {
    if(ex.isCausedByNetworkIssue()) {
      // 何らかの事情でネットワークに接続できませんでした
    } else {
      int statusCode = ex.getStatusCode();
      switch(statusCode) {
        case 400:  // Bad Request
          // リクエストに不正があるかレートリミットに達しました
          break;
        case 401:  // Unauthorized
          // 認証情報に間違いがあります
          break;
        case 403:  // Forbidden
          // リクエストが拒否されました
          break;
        case 404:  // Not Found
          // 存在しないリソースにアクセスしました
          break;
        case 500:  // Server Internal Error
          // サーバ障害が発生しています
          break;
        case 502:  // Bad Gateway
          // Twitterサービスに障害が発生しています
          break;
        case 503:  // Service Unavailable
          // Twitterが過負荷状態です
          break;
      }
    }
  }
  
  noLoop();
}

void draw() {
  // 
}

実行すると、コンソールにフレンド一覧が列挙されます。

ただこれ、けっこう時間がかかります。ですので、フォロー一覧の読み込み処理は別スレッドに移した方がよいでしょう。適当にマルチスレッド化するとこうなります（本当に適当なので、よい子はマネしないでね）。

【ソースコード（展開してご覧ください）】

import twitter4j.*;

final String USER_NAME = "tercel_s";

Twitter twitter;
boolean isLoading;

void setup() {
  twitter = new TwitterFactory().getInstance();
  
  // マルチスレッドだよ∩( ・ω・)∩
  new Thread() {
    public synchronized void run() {
      isLoading = true;
      try {
        // 指定したユーザがフォローしているアカウント一覧を取得
        IDs friends = twitter.getFriendsIDs(USER_NAME, -1);
        long[] friendsIDs = friends.getIDs();
        
        // Twitter ID からユーザ情報を逐次取得して表示
        for(long id : friendsIDs) {
          User user = twitter.showUser(id);
          println(user.getName() + " (" + user.getScreenName() + ")");
        }
        isLoading = false;
      } catch (TwitterException ex) {
        /* 省略 */
      } finally {
        isLoading = false;  // デッドロック回避用
      }
    }
  }.start();
}

void draw() { 
  // 
}

これさえできれば、デモ動画の作品を作るのは簡単です。ほとんど昨日のパーティクルをユーザのアイコンに差し換えただけです。

【おまけソースコード（展開してご覧ください）】

import processing.opengl.*;
import javax.media.opengl.*;
import twitter4j.*;

final int FIELD_SIZE = 1000;
final int FIELD_STEP =  100;

final String USER_NAME = "tercel_s";  
boolean      isLoading;

List<Agent> agentList;  // フレンドアイコンのリスト

long time;

void setup() {
  size(800, 600, OPENGL);  
  
  agentList = new ArrayList<Agent>();

  // マルチスレッドでフレンドを読み込む
  new Thread(new FriendLoader()).start();
  hint(ENABLE_DEPTH_SORT);
}

void draw() { 
  background(0);
  
  /* ================== */
  /* カメラを適当に設定 */
  /* ================== */
  float angle = 0.05f * radians(time);
  
  camera(800, -200,   0, 
           0,    0,   0, 
           0,    1,   0);
  rotateY(angle);  // ちょっと回す

  pushStyle();
  
  // テクスチャの加算合成を有効にする
  PGraphicsOpenGL pgl = (PGraphicsOpenGL)g;
  GL gl = pgl.beginGL();
  
  gl.glEnable(GL.GL_BLEND);
  gl.glBlendFunc(GL.GL_SRC_ALPHA, GL.GL_ONE_MINUS_SRC_ALPHA);
  pgl.endGL();  
  
  noStroke();
  fill(0);
  
  // フレンド一覧を表示
  for(int i = 0; i < agentList.size(); ++i) {
    agentList.get(i).update();
  }

  popStyle();

  /* ============== */
  /* 地面とかの描画 */
  /* ============== */
  pushStyle();
  // 地面の描画 
  stroke(0, 255, 0, 100);
  for(int i = -FIELD_SIZE; i <= FIELD_SIZE; i += FIELD_STEP) { 
    line(          i, 0, -FIELD_SIZE,          i, 0, FIELD_SIZE);   
    line(-FIELD_SIZE, 0,           i, FIELD_SIZE, 0,          i);
  }
  
  popStyle();
  
  time++;
  mm.addFrame();
  if(time > 60 * 60) {
    mm.finish();
    background(0);
    noLoop();
  }
}


// フレンドを読み込むためのクラス
class FriendLoader implements Runnable {
  public void run() {
    isLoading = true;
    Twitter twitter = new TwitterFactory().getInstance();
    try {
      // まずは自分自身をリストに追加
      User me = twitter.showUser(USER_NAME);
      PImage myIcon = loadImage(me.getProfileImageURL().toString());
      myIcon.resize(100, 100);
      agentList.add(new Agent(myIcon, new PVector()));
            
      // 自分がフォローしているアカウント一覧を取得
      IDs friends = twitter.getFriendsIDs(USER_NAME, -1);
      long[] friendsIDs = friends.getIDs();
        
      // Twitter ID からユーザ情報を逐次取得
      for(long id : friendsIDs) {
        User user = twitter.showUser(id);
          
        // アイコン取得とリサイズ
        PImage profileImage = loadImage(user.getProfileImageURL().toString());
        if(profileImage == null) continue;
        
        profileImage.resize(100, 100);

        // 位置の設定
        float x = random(-1000, 1000);
        float z = random(-1000, 1000);
        agentList.add(new Agent(profileImage, new PVector(x, 0, z)));
      }
      
      isLoading = false;
    } catch (TwitterException ex) {
      /* 省略 */
    } finally {
      isLoading = false;  // デッドロック回避用
    }
  }
}


// エージェント
class Agent {
  private final float TERRITORY_SIZE = 2000;
  private final float NOISE_SCALE    = 0.001f;
  private final float AGENT_SIZE     = 50; 
  
  private PImage  tex;     // テクスチャ
  
  private PVector center;  // 縄張りの中心座標
  private PVector offset;  // 中心からのオフセット
  private PVector pos;     // 座標
  
  private long time;      // 経過時間
  float xOffset, zOffset; // ノイズ生成用
  
  Agent(PImage tex, PVector center) {
    this.center = center;
    this.tex    = tex;
    pos         = new PVector();
    offset      = new PVector();
    xOffset     = random(TERRITORY_SIZE);
    zOffset     = random(TERRITORY_SIZE);
  }
  
  PVector getPos() {
    return pos;
  }
  
  void update() {
    if (tex == null) return;
    time++;
    
    offset.x = (noise((time + xOffset) * NOISE_SCALE) - 0.5f) * TERRITORY_SIZE;
    offset.z = (noise((time + zOffset) * NOISE_SCALE) - 0.5f) * TERRITORY_SIZE;
    
    pushMatrix();
    pos.set(center.x + offset.x, center.y, center.z + offset.z);
    translate(pos.x, pos.y, pos.z);
    
    // モデル-ビュー行列を取得
    PMatrix3D builboardMat = (PMatrix3D)g.getMatrix();
    // 回転成分を単位行列に
    builboardMat.m00 = builboardMat.m11 = builboardMat.m22 = 1;
    builboardMat.m01 = builboardMat.m02 = 
    builboardMat.m10 = builboardMat.m12 = 
    builboardMat.m20 = builboardMat.m21 = 0;
  
    resetMatrix();
    applyMatrix(builboardMat);
        
    beginShape(QUADS);
    texture(tex);
    vertex(-0.5f * AGENT_SIZE, -AGENT_SIZE, 0, 0,         0);
    vertex(-0.5f * AGENT_SIZE,           0, 0, 0,         tex.height-1);
    vertex( 0.5f * AGENT_SIZE,           0, 0, tex.width, tex.height-1);
    vertex( 0.5f * AGENT_SIZE, -AGENT_SIZE, 0, tex.width, 0);    
    endShape();
        
    popMatrix();
  }
}

あとはこの 2 つを適当に組み合わせれば、冒頭の動画のようなものが作れると思います。

なしくずしてきにめでたしめでたし。

そんなこんなで久しぶりに Twitter ネタを書きましたが、このプログラムにはまだまだアラがあります。

まず、Twitter の API 呼び出しには制限が設けられており、このプログラムを試しているとあっという間に API 上限に引っかかってしまいます。フォローしているアカウントの情報を一人ずつ取得していくと、アホみたいなスピードで限界に近付いていきます。

次に、計算機資源の問題があります。今はとりあえずフレンド情報を全読みしていますが、それこそ何十万人もフォローしているような人が動かしたら、あっという間にメモリを喰い尽くして死ぬでしょう。

それでもまぁ、とりあえず Twitter と P5 をくっつけてそれっぽい事ができたので今日はよしとします。ふぅ、つかれた(´-ω-`)

Processingできらきら3次元パーティクル

2012-02-14T20:37:00.003+09:00

今日はこういうものを作ります。たぶんきれいですね。

このスケッチは、いくつかのテクニックを組み合わせて作ります。

恐らく、Processing ユーザにとっては馴染みの薄いものもあるかと思いますので、今日はそれらをひとつひとつ紹介しつつ完成形に持っていこうと考えています。

【ビルボードをつくる】

ビルボードとは、3次元空間中のカメラに対して常に正面を向くような2次元オブジェクトの事です。一昔前のゲームなどではよく応用されていた技術です（今もかな）。これ自体が独立したひとつのトピックなので、まずはここから始めてみましょう。

以下のサンプルは、テキスト（Hello, World）をビルボード上に表示するというものです。

【理論の解説】

ビルボードは、以下のような同次形のモデル・ビュー変換行列
\begin{align}
{}^{\Sigma_{\mathrm{View}}}\vect{M}_{\Sigma_{\mathrm{Model}_{i}}} =
\left(\!
\begin{array}{ccc|c}
m_{00} & m_{01} & m_{02} & m_{03} \\
m_{10} & m_{11} & m_{12} & m_{13} \\
m_{20} & m_{21} & m_{22} & m_{23} \\ \hline
m_{30} & m_{31} & m_{32} & m_{33}
\end{array}
\!\right)
\end{align}――のうち、左上の3×3の小行列を単位行列にした結果
\begin{align}
\left( \!
\begin{array}{ccc|c}
1 & 0 & 0 & m_{03} \\
0 & 1 & 0 & m_{13} \\
0 & 0 & 1 & m_{23} \\ \hline
m_{30} & m_{31} & m_{32} & m_{33}
\end{array}
\!\right)
\end{align}を、新たなモデル・ビュー行列として設定します。たったこれだけ。かんたんですね。

【ソースコード（展開してご覧ください）】

import processing.opengl.*;

PFont font;
long  time;

void setup() {
  size(800, 600, OPENGL);
  
  // フォント生成（メイリオ）
  font = createFont("Meiryo", 48);
  
  // デプスソートを有効にする
  // これがないとテキスト表示がおかしくなる
  hint(ENABLE_DEPTH_SORT);
  
}

void draw() {
  background(0);
  
  /* ================== */
  /* カメラを適当に設定 */
  /* ================== */
  float angle = 0.5f * radians(time);
  
  camera(0, -300*sin(angle), 300, 
         0,               0,   0, 
         0,               1,   0);
  rotateY(angle);  // ちょっと回す

  /* ================ */
  /* ビルボーディング */
  /* ================ */
  pushMatrix();
  translate(100, 0, 0);  // 物体を並進（ローカル座標系）
  
  // モデル-ビュー行列を取得
  PMatrix3D builboardMat = (PMatrix3D)g.getMatrix();
  // 回転成分を単位行列に
  builboardMat.m00 = builboardMat.m11 = builboardMat.m22 = 1;
  builboardMat.m01 = builboardMat.m02 = 
  builboardMat.m10 = builboardMat.m12 = 
  builboardMat.m20 = builboardMat.m21 = 0;

  resetMatrix();
  applyMatrix(builboardMat);
  
  fill(255);
  stroke(0);
  textFont(font, 24);
  textAlign(CENTER);
  text("Hello, World!", 0, 0);

  popMatrix();

  /* ============== */
  /* 地面とかの描画 */
  /* ============== */
  pushStyle();
  // 地面の描画 
  stroke(0, 255, 0);
  for(int x = -1000; x <= 1000; x += 100)
    line(x, 0, -1000, x, 0, 1000); 
  for(int z = -1000; z <= 1000; z += 100)
    line(-1000, 0, z, 1000, 0, z);
  
  // 回転軸の描画
  stroke(255, 0, 0);
  line(0, -500, 0, 0, 500, 0);
  popStyle();
  
  time++;
}

【自然なランダムウォークをつくる】

自然なランダムウォークというと少しおかしな日本語ですが、あちこちに動きが飛ばず、なめらかな軌跡を描くような動きを作ります。

タネを明かせば Perlin ノイズを何のひねりもなく使っただけで、当たり判定とかは考慮していません。

【ソースコード（展開してご覧ください）】

ここでは、動く個体を Agent クラスにまとめています。

import processing.opengl.*;

List<Agent> agentList;
long  time;

void setup() {
  size(800, 600, OPENGL);
  
  agentList = new ArrayList<Agent>();
  for(int i = 0; i < 200; ++i) {
    float x = random(-200, 200);
    float z = random(-200, 200);
    agentList.add(new Agent(new PVector(x, 0, z)));
  }
}

void draw() {
  background(0);
  camera();
  lights();
    
  /* ================== */
  /* カメラを適当に設定 */
  /* ================== */
  float angle = 0.1f * radians(time);
  
  camera(800, -200,   0, 
         0,      0,   0, 
         0,      1,   0);
  rotateY(angle);  // ちょっと回す

  pushStyle();
  noStroke();
  colorMode(HSB, agentList.size() -1, 100, 100);
  for(int i = 0; i < agentList.size(); ++i) {
    fill(i, 100, 100);
    agentList.get(i).update();
  }  
  popStyle();

  /* ============== */
  /* 地面とかの描画 */
  /* ============== */
  pushStyle();
  // 地面の描画 
  stroke(0, 255, 0);
  for(int x = -1000; x <= 1000; x += 100)
    line(x, 0, -1000, x, 0, 1000); 
  for(int z = -1000; z <= 1000; z += 100)
    line(-1000, 0, z, 1000, 0, z);
  
  // 回転軸の描画
  stroke(255, 0, 0);
  line(0, -500, 0, 0, 500, 0);
  popStyle();
  
  time++;
}

class Agent {
  private final float TERRITORY_SIZE = 2000;
  private final float NOISE_SCALE    = 0.001f;
  private final float AGENT_SIZE     = 10; 
  
  private PVector center;  // 縄張りの中心座標
  private PVector offset;  // 中心からのオフセット
  private PVector pos;     // 座標
  
  private long time;      // 経過時間
  float xOffset, zOffset; // ノイズ生成用
  
  Agent(PVector center) {
    this.center = center;
    pos         = new PVector();
    offset      = new PVector();
    xOffset     = random(TERRITORY_SIZE);
    zOffset     = random(TERRITORY_SIZE);
  }
  
  PVector getPos() {
    return pos;
  }
  
  void update() {
    time++;
    
    offset.x = (noise((time + xOffset) * NOISE_SCALE) - 0.5f) * TERRITORY_SIZE;
    offset.z = (noise((time + zOffset) * NOISE_SCALE) - 0.5f) * TERRITORY_SIZE;
    
    pushMatrix();
    pos.set(center.x + offset.x, center.y - AGENT_SIZE * 0.5, center.z + offset.z);
    translate(pos.x, pos.y, pos.z);
    
    box(AGENT_SIZE);
        
    popMatrix();
  }
}

【合体させる】

あとは、この2つのテクニックをてきとうに合体させるだけです。

下図（左）のように、とりあえずランダムに動き回るビルボードを山のように作って、それにテクスチャを貼れば完成。意外と簡単ですね。

左：ビルボードの輪郭 / 右：テクスチャを貼った結果

ここで、glBlendFunc メソッドを使って加算合成モードにすると、ビルボード同士が重なった際に輝度が上がります。完成形のソースコードを以下に示します。

【ソースコード（展開してご覧ください）】

import processing.opengl.*;
import javax.media.opengl.*;

PGraphics   tex;
List<Agent> agentList;
long        time;

void setup() {
  size(800, 600, OPENGL);
 
  // テクスチャ生成
  tex = createTexture();
  
  agentList = new ArrayList<Agent>();
  for(int i = 0; i < 1000; ++i) {
    float x = random(-1000, 1000);
    float z = random(-1000, 1000);
    agentList.add(new Agent(new PVector(x, 0, z)));
  }
}

void draw() {
  background(0);
    
  /* ================== */
  /* カメラを適当に設定 */
  /* ================== */
  float angle = 0.05f * radians(time);
  
  camera(800, -200,   0, 
           0,    0,   0, 
           0,    1,   0);
  rotateY(angle);  // ちょっと回す

  pushStyle();
  
  // テクスチャの加算合成を有効にする
  PGraphicsOpenGL pgl = (PGraphicsOpenGL)g;
  GL gl = pgl.beginGL();
  gl.glDisable(GL.GL_DEPTH_TEST);
  gl.glEnable(GL.GL_BLEND);
  gl.glBlendFunc(GL.GL_SRC_ALPHA, GL.GL_ONE);
  pgl.endGL();  
  
  noStroke();
  fill(0);
  for(int i = 0; i < agentList.size(); ++i) {
    agentList.get(i).update();
  }

  popStyle();

  /* ============== */
  /* 地面とかの描画 */
  /* ============== */
  pushStyle();
  // 地面の描画 
  stroke(0, 255, 0, 100);
  for(int x = -1000; x <= 1000; x += 100)
    line(x, 0, -1000, x, 0, 1000); 
  for(int z = -1000; z <= 1000; z += 100)
    line(-1000, 0, z, 1000, 0, z);
  
  popStyle();
  
  time++;
}

PGraphics createTexture() {
  PGraphics tex = createGraphics(100, 100, P2D);
  float r = tex.width * 2 / 3;
  
  tex.beginDraw();
  tex.background(0);
  tex.noStroke();
  tex.fill(100, 150, 255);
  tex.ellipse(tex.width/2, tex.height/2, r, r);
  tex.filter(BLUR, 10.0);
  tex.endDraw();
  return tex;
}

class Agent {
  private final float TERRITORY_SIZE = 2000;
  private final float NOISE_SCALE    = 0.001f;
  private final float AGENT_SIZE     = random(50); 
  
  private PVector center;  // 縄張りの中心座標
  private PVector offset;  // 中心からのオフセット
  private PVector pos;     // 座標
  
  private long time;      // 経過時間
  float xOffset, zOffset; // ノイズ生成用
  
  Agent(PVector center) {
    this.center = center;
    pos         = new PVector();
    offset      = new PVector();
    xOffset     = random(TERRITORY_SIZE);
    zOffset     = random(TERRITORY_SIZE);
  }
  
  PVector getPos() {
    return pos;
  }
  
  void update() {
    time++;
    
    offset.x = (noise((time + xOffset) * NOISE_SCALE) - 0.5f) * TERRITORY_SIZE;
    offset.z = (noise((time + zOffset) * NOISE_SCALE) - 0.5f) * TERRITORY_SIZE;
    
    pushMatrix();
    pos.set(center.x + offset.x, center.y, center.z + offset.z);
    translate(pos.x, pos.y, pos.z);
    
    // モデル-ビュー行列を取得
    PMatrix3D builboardMat = (PMatrix3D)g.getMatrix();
    // 回転成分を単位行列に
    builboardMat.m00 = builboardMat.m11 = builboardMat.m22 = 1;
    builboardMat.m01 = builboardMat.m02 = 
    builboardMat.m10 = builboardMat.m12 = 
    builboardMat.m20 = builboardMat.m21 = 0;
  
    resetMatrix();
    applyMatrix(builboardMat);
        
    beginShape(QUADS);
    texture(tex);
    vertex(-0.5f * AGENT_SIZE, -AGENT_SIZE, 0, 0,         0);
    vertex(-0.5f * AGENT_SIZE,           0, 0, 0,         tex.height-1);
    vertex( 0.5f * AGENT_SIZE,           0, 0, tex.width, tex.height-1);
    vertex( 0.5f * AGENT_SIZE, -AGENT_SIZE, 0, tex.width, 0);    
    endShape();
        
    popMatrix();
  }
}

【予告】 Processingで講演します（？）

2012-02-13T13:42:00.003+09:00

少し先のお話ですが、来たる3月20日～23日に開催される電子情報通信学会・総合大会に参戦する事になりました。講演分野は『パターン認識・メディア理解』だそうです。

今回は「Processing で『立体おえかきツール』を作ってみたよ」みたいな内容です。とりあえず予告篇のようなデモムービーをてきとうに作ったので、こっそりご紹介しておきますね（※注意：音が鳴ります）。

あまり大規模な制作物ではありませんが、もしよろしければどうぞおたのしみに。

【関連研究】

Robert C. Zeleznik, Kenneth P. Herndon, and John F. Hughes, “SKETCH: An interface for sketching 3D scenes”, ACM Transactions on Graphics (SIGGRAPH Proceedings), pp.163-170, 1996.
Takeo Igarashi, Satoshi Matsuoka, Hidehiko Tanaka, “Teddy: A Sketching Interface for 3D Freeform Design”, ACM Transactions on Graphics (SIGGRAPH Proceedings), pp.409-416, 1999.
Dorsey Julie, Xu Songhua, Smedresman Gabe, Rushmeir Holly, McMillan Leonard, “The Mental Canvas: A Tool for Conceptual Architectural Design and Analysis”, Proceedings of Pacific Graphics, pp.201-210, 2007.
Schmidt Ryan, Khan Azam, Singh Karan, Kurtenbach Gord, “Analytic Drawing of 3D Scaffolds”, ACM SIGGRAPH Asia 2009 papers 28(5), pp.149:1-149:10, 2009.
“SKETCHUP”, Google, 2009.

……などなど。

Arduinoはじめます

2012-02-09T22:23:00.001+09:00

研究が一段落して論文を書き終えてしまったためか、なんとも言えない虚脱感に襲われています。

そこで、新たな生きがいを探すべく、Arduino というおもちゃを買ってきました。

入門用の『Arduinoをはじめようキット』

Arduino は、自分で回路を組み、マイコンにプログラムを送り込むと、とってもいいものができる素敵な装置なのだそうです（※ 厳密には装置の名称ではないのだけど）。

電気電子な人じゃなくても電子工作ごっこができ、マイコンボード本体も中高生のお小遣いで買えてしまう程のお手頃価格。これなら僕にもできるかな。

【開封の儀式】

ちなみに内容物は以下の通りでした。

Arduino ボード（Arduino UNO）
ブレッドボード
ジャンパワイヤ
モメンタリ型プッシュボタンスイッチ（タクトスイッチ）
CdS セル（光センサ）
抵抗（10kΩ×2、1kΩ×3、330Ω×1）

電気回路や電子回路は専門外すぎてほとんど分かりません。ちょっとずつ慣れていこう……。

【タクトスイッチの実験】

とりあえず基本という事で、スイッチを使って LED を点けたり消したりできるように回路を作ってみました。

押しているときだけ光ります。

ちなみに全く同一の回路でも、プログラムを書き換えると挙動を変える事が可能です。たとえば、スイッチから指を離しても LED の点灯状態が維持されるようにできます。

【CdS セルの実験】

続いて、タクトスイッチを CdS セルに接ぎ換えてみる事に。これは、光に反応して抵抗値が変動するおもしろいセンサだそうです。

実は最初、どれだけ CdS セルに光を当てても何故か LED が点灯せず、『初期不良だな、ぷんぷん』と思ったら、本来 1kΩ の抵抗を使うべきところに 10kΩ の抵抗を使っていた事が発覚(・ω・lll)

とんでもない凡ミスでした。

【3色 LED の実験】

最後に、3色の LED を使ってグラデーションを作ろうとして失敗してみました。

HSV 表色系の H をじわじわ変動させたらグラデーションっぽい事ができるんじゃね？と思ったのですが、なんか残念なものができただけでした。

……というわけで、一応セットされている部品はひととおり使ってみましたが、これらを組み合わせて何ができるのかは全くの未知数です。

慣れてきたら Processing とも組み合わせてみたいなぁ。

※ 今日の日記に使用した写真・動画はデジタルハリネズミで撮影しました。なんかレトロっぽいよね。

C++と仲直り計画： Cの構造体をスマートポインタで扱う（同時上映：ラムダ式もほんのちょっとかじる）

2012-01-28T18:20:00.000+09:00

いきなりですが今日は boost::shared_ptr を構造体に対して用いるお話です。ここまで読んで、『何を今さら』と思った方はブラウザをそっと閉じると幸せになれますよ。

あと、今日の日記は解説というよりは備忘録的な色が強いので、僕と波長が合わないとなかなか内容に共感して頂けないかも知れません。そういう方もブラウザをそっｔ（ｒｙ

【だらだら C 言語からスマート C++ へ】

僕が C / C++ に苦手意識を抱く一つの理由がメモリ管理です。まず、動的に確保したメモリをどのタイミングで解放するか、あるいは解放のし忘れはないかを考えるのが煩わしい。次に、メモリ解放用のコードを書くのが面倒。

百聞は一見に何とやらで、メモリを動的に確保 / 解放するコードの例を以下に示します。
（ちなみに OpenCV を使っているのは個人的な趣味です）（ごめんね）

【OpenCV で画像領域を動的確保/解放（生ポインタ）】

#include<iostream>

#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>

#ifdef _DEBUG
#pragma comment(lib, "opencv_core231d.lib")
#pragma comment(lib, "opencv_highgui231d.lib")
#else
#pragma comment(lib, "opencv_core231.lib")
#pragma comment(lib, "opencv_highgui231.lib")
#endif

using namespace std;

int main() {
    // 400 × 300 で画像を作成
    IplImage* img = cvCreateImage(cvSize(400, 300), IPL_DEPTH_8U, 3);

    // img を表示
    cvShowImage("作った画像を見せびらかすウィンドウ", img);
    cv::waitKey(0);

    // img を解放
    cvReleaseImage(&img);
    img = NULL;
    
    cout << "メモリを適当に解放しました\n";
    return 0;
}

念のため、上のコードを簡単に説明すると、cvCreateImage で動的にメモリを確保して、cvReleaseImage で解放しています。これ、IplImage 型の変数が 1 つだけしかないからよいものの、変数の数だけ解放用コードが並ぶのはさぞ目障りなものでしょう。

ところが、boost というライブラリには、動的に生成されたオブジェクトを自動的にデストラクトしてくれる便利なスマートポインタが用意されているのです。

スマートポインタにもいくつか種類がありますが、今回は最も汎用性が高い（と個人的に考えている）boost::shared_ptr を使ってみましょう。

まず、boost::shared_ptr<iplimage> 用のカスタムデリータをファンクタの形で定義します（ここでは ImageReleaser と名付けた）。やっている事は、単に関数呼び出し演算子をオーバーロードして、cvReleaseImage を実行しているだけです。簡単ですね。

次に、このデリータオブジェクトを boost::shared_ptr<iplimage> のコンストラクタ（第2引数）に渡してやります。すると、たったこれだけで生成されたオブジェクトが自動解放されるようになります。やったね！

【OpenCV で画像領域を動的確保/解放（スマートポインタ）】

#include<iostream>

#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>

#include<boost/shared_ptr.hpp>

#ifdef _DEBUG
#pragma comment(lib, "opencv_core231d.lib")
#pragma comment(lib, "opencv_highgui231d.lib")
#else
#pragma comment(lib, "opencv_core231.lib")
#pragma comment(lib, "opencv_highgui231.lib")
#endif

using namespace std;

// 解放用のファンクタ
class ImageReleaser {
public:
    void operator()(IplImage* img) {
        cvReleaseImage(&img);
        cout << "メモリを適当に解放しました\n";
    }
};


int main() {
    
    typedef boost::shared_ptr<IplImage> ImgPtr;
    
    // 400 × 300 で画像を作成
    ImgPtr img(cvCreateImage(cvSize(400, 300), IPL_DEPTH_8U, 3), 
               ImageReleaser());    // ←解放用ファンクタを指定

    // img を表示
    cvShowImage("作った画像を見せびらかすウィンドウ", img.get());
    cv::waitKey(0);

    // img は自動的に解放される
    return 0;
}

おー、できた。でもまだちょっとアレですね。

プログラム中で一度しか使われないデリータごときのためにわざわざファンクタを定義するのは大仰です。それに、極めて限定的な用途の小さなクラスがソース内に散らかるのはどうも汚く思えます。

ここで、先ほどの実装の改善すべき点を洗ってみましょう。

デリータは基本的に一度きりの使い捨てなので名前など要らないでしょうし、必要な時・必要な場所に実装を書けば「小さなクラスがソース内に散らかる」といった気持ち悪さもなくなります。

これを実現するのが boost::lambda です。そう、みんな大好きラムダ式のライブラリですね。ここではあまりマニアックな使い方はせず、単に邪魔なファンクタをこの世から消し去る目的のためだけに導入しました。

先ほどのソースを書き換えると以下のようになります。

【OpenCV で画像領域を動的確保/解放（スマートポインタ + ラムダ式）】

#include<iostream>

#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>

#include<boost/shared_ptr.hpp>

#include<boost/lambda/bind.hpp>
#include<boost/lambda/lambda.hpp>

#ifdef _DEBUG
#pragma comment(lib, "opencv_core231d.lib")
#pragma comment(lib, "opencv_highgui231d.lib")
#else
#pragma comment(lib, "opencv_core231.lib")
#pragma comment(lib, "opencv_highgui231.lib")
#endif

using namespace std;

int main() {
    using namespace boost::lambda;
    typedef boost::shared_ptr<IplImage> ImgPtr;
    
    // 400 × 300 で画像を作成
    ImgPtr img(cvCreateImage(cvSize(400, 300), IPL_DEPTH_8U, 3),
        (bind(cvReleaseImage, &_1),
            cout << constant("メモリを適当に解放しました\n")));

    // img を表示
    cvShowImage("作った画像を見せびらかすウィンドウ", img.get());
    cv::waitKey(0);

    // img は自動的に解放される
    return 0;
}

この例のほかにも、ラムダ式を使えばソート関数にカスタムコンパレータを与えるコードをスマートに書けるでしょう（3次元 Delaunay 分割では boost レスで作ったので、わざわざコンパレータをファンクタにしていた）。

ちなみに OpenCV には cv::Ptr というすごい便利なスマートポインタが用意されています。というか、僕は boost を入れる前にもっぱらこちらを使っていました。

これは IplImage などのカスタムデリータが自動的に適用されるので、それを使えばひとまず上記の議論はことごとく不要になるのですが、たとえば悪名高い COM オブジェクトとか、C ベースで書かれた過去の重要な資産とかを安全に扱いたい場合には知っておくと幸せになれるような気がします。

Processing.jsで作った「ひこうき」を高速化してみた

2012-01-25T16:40:00.000+09:00

以前、Processing.js で『3D ひこうき』を作りました。

これは P3D モードのサポートがアレな Processing.js 環境でも、なんとか 3D っぽいグラフィックを出そうというチャレンジだったのですが、速度的な意味で非常に残念なものができあがっただけでした。

今回は、ループの中でオブジェクトを new しないように心がけたり、描画関数の呼び出しを削減したりといろいろ高速化を試みました（→ OpenProcessing）。

うーん……。

僕の iPhone 3GS だとまだカクカクしています。

もうちょっと形状を単純にするとかなり速く動く気がするのですが……。

Processingでてきとう3次元バンプマッピング

2012-01-23T17:30:00.000+09:00

学位論文の方が一段落しそうなので、空き時間を使ってバンプマッピング（のようなもの）を実装してみました。ちょっとしたリハビリです。

バンプマッピングとは、ポリゴン上の凹凸を擬似的に再現する技術で、テクスチャマッピングでは表現できない質感表現を実現できます。

実際に動作するデモは OpenProcessing で公開していますので、是非ご覧ください。

三次元・直線・最小自乗法（vs. ぷで）

2012-01-19T16:04:00.005+09:00

ぷでさんが面白い記事を書いていたので、僕も同じテーマで書いてみる事にしたよ。

「3次元における直線 \(L_{i}\: \left( i = 0,1,\cdots, n-1 \right) \) の尤もらしい交点 \( \vect{p} \) を最小自乗法で求める」

これは、位置と姿勢が既知である多視点画像における対応点3次元位置を推定する際にしばしば用いられる議論だ（下図参照）。3次元空間中では一般的に複数の直線の交点は一意には定まらず、いわゆる “ねじれ” の関係になる。その際、できるだけ誤差を抑えるために、各直線が最接近する点を採用する事が望ましい。

複数カメラの視線を基に被写体の3次元座標を得るイメージ．
カメラ光学中心から各対応点を貫く直線の交点が3次元位置となる．

直線 \(L_{i}\) が通る点を \( \vect{b}_{i} \) 、傾き（方向）を \( \vect{d}_{i} \) とすると、\(L_{i}\) 上の任意の点はパラメータ \(t_{i}\) を用いて
\begin{align}
L_{i}(t_{i}) = t_{i}\vect{d}_{i} + \vect{b}_{i}
\end{align}と表す事ができる。

また、直線上の座標 \(L_{i}(t_{i})\) と点 \(\vect{p}\) の距離の自乗 \(f_{i}(t_{i},\,\vect{p})\) は以下のようになる。
\begin{align}
f_{i}(t_{i},\,\vect{p}) &= \left| L_{i}(t_{i}) - \vect{p} \right|^{2} \\
&= \left| t_{i}\vect{d}_{i} + \vect{b}_{i} - \vect{p} \right|^{2}
\end{align}

ここで、すべての直線上の点 \(L_{i}(t_{i})\) から \(\vect{p}\) までの距離自乗の総和を
\begin{align}
& F(t_{0},\, \cdots,\, t_{i-1},\, t_{i},\, t_{i+1},\, \cdots, \, t_{n-1}, \, \vect{p}) \\
=& \frac{1}{2} \sum_{i=0}^{n-1} \left| t_{i} \vect{d}_{i} + \vect{b}_{i} - \vect{p} \right|^{2} \\
=& \frac{1}{2} \sum_{i=0}^{n-1} \left( t_{i} \vect{d}_{i} + \vect{b}_{i} - \vect{p} \right)^{t} \left( t_{i} \vect{d}_{i} + \vect{b}_{i} - \vect{p} \right)
\end{align}と置き、関数 \(F(t_{0},\, \cdots,\, t_{i-1},\, t_{i},\, t_{i+1},\, \cdots, \, t_{n-1}, \, \vect{p})\) の最小化問題を解く。

まず、上式を各 \( t_{i} \) で偏微分して \( 0 \) と置き、
\begin{align}
& \frac{\partial }{\partial t_{i}} F(t_{0},\, \cdots,\, t_{i-1},\, t_{i},\, t_{i+1},\, \cdots, \, t_{n-1}, \, \vect{p}) \\
=& \frac{\partial}{\partial t_{i}} \left[ \frac{1}{2} \sum_{i=0}^{n-1} \left( t_{i} \vect{d}_{i} + \left( \vect{b}_{i} - \vect{p} \right) \right)^{t} \left( t_{i} \vect{d}_{i} + \left( \vect{b}_{i} - \vect{p} \right) \right) \right] \\
=& \frac{\partial}{\partial t_{i}} \left[ \frac{1}{2} \sum_{i=0}^{n-1} t_{i}^{2}\vect{d}_{i}^{t}\vect{d}_{i} + \sum_{i=0}^{n-1} t_{i} \vect{d}_{i}^{t} \left( \vect{b}_{i} - \vect{p} \right) + \frac{1}{2} \sum_{i=0}^{n-1} \left( \vect{b}_{i} - \vect{p} \right)^{t} \left( \vect{b}_{i} - \vect{p} \right) \right] \\
=& t_{i}\vect{d}_{i}^{t}\vect{d}_{i} + \vect{d}_{i}^{t} \left( \vect{b}_{i} - \vect{p} \right) = 0
\end{align}
次に \( \vect{p} \) で偏微分して同様に \( 0 \) と置くと
\begin{align}
&\frac{\partial }{\partial \vect{p}} F(t_{0},\, \cdots,\, t_{i-1},\, t_{i},\, t_{i+1},\, \cdots, \, t_{n-1}, \, \vect{p}) \\
=& \frac{\partial}{\partial \vect{p}} \left[ \frac{1}{2} \sum_{i=0}^{n-1} \left( \left( t_{i} \vect{d}_{i} + \vect{b}_{i} \right) - \vect{p} \right)^{t} \left( \left( t_{i} \vect{d}_{i} + \vect{b}_{i} \right) - \vect{p} \right) \right] \\
=& \frac{\partial}{\partial \vect{p}} \left[ \frac{1}{2} \sum_{i=0}^{n-1} \left( t_{i} \vect{d}_{i} + \vect{b}_{i} \right)^{t} \left( t_{i} \vect{d}_{i} + \vect{b}_{i} \right) - \sum_{i=0}^{n-1} \left( t_{i} \vect{d}_{i} + \vect{b}_{i} \right)^{t} \vect{p} + \frac{1}{2}\sum_{i=0}^{n-1} \vect{p}^{t}\vect{p} \right] \\
=& -\sum_{i=0}^{n-1} \left( t_{i} \vect{d}_{i} + \vect{b}_{i} \right) + n\vect{p} = \vect{0}
\end{align}が得られる。

これらをまとめると、未知数 \( t_{0},\, \cdots,\, t_{n-1},\, \vect{p} \) に関する次の方程式になる。
\begin{align}
\begin{pmatrix}
\vect{d}_{0}^{t}\vect{d}_{0} & 0 & \cdots & 0 & -\vect{d}_{0}^{t} \\
0 & \vect{d}_{1}^{t}\vect{d}_{1} & \ddots & \vdots & -\vect{d}_{1}^{t} \\
\vdots & \ddots & \ddots & 0 & \vdots \\
0 & \cdots & 0 & \vect{d}_{n-1}^{t}\vect{d}_{n-1} & -\vect{d}_{n-1}^{t} \\
-\vect{d}_{0} & -\vect{d}_{1} & \cdots & -\vect{d}_{n-1} & \vect{E}_{3,\,3}
\end{pmatrix}
\begin{pmatrix}
t_{0} \\
t_{1} \\
\vdots \\
t_{n-1} \\
\vect{p}
\end{pmatrix} =
\begin{pmatrix}
-\vect{d}_{0}^{t} \vect{b}_{0} \\
-\vect{d}_{1}^{t} \vect{b}_{1} \\
\vdots \\
-\vect{d}_{n-1}^{t} \vect{b}_{n-1} \\
\sum_{i=0}^{n-1}\vect{b}_{i}
\end{pmatrix}
\end{align}

ここで、\( \vect{E}_{3,\,3} \) は 3 階の単位行列である。上記の方程式を LU 分解か何かを使って解けばたぶん最急降下法を用いなくても \( \vect{p} \) を得る事ができる（…んじゃないかな）。

ちなみに boost を用いた LU 分解の実装例はこのへん。OpenCV による実装例はこのへん。

素晴らしき日々へ

2012-01-08T21:28:00.000+09:00

カルボナーラを食べるとカルボニル基という言葉を思い出す。

ボロネーゼを食べると、ドロネー図を思い出す（ついでにボロノイ図も思い出す）。

普段からこのように突飛な方向へと思考回路がスイッチされてしまうため、話し相手はそれなりに苦慮しているらしい。申し訳ない。

今日は親友の祥月命日だ。

本当は一緒に連れて逝って欲しかったのだが、僕は生きることになった。

死ぬというのはどういう事なのだろう。僕はあまり宗教への頓着がないせいか、恐らく『生まれる前と同じ世界』に還るのだろうと考えている。

最近の気分転換は専らマリオだ。

小さい頃から高いところがダメで、たとえゲームであっても足が竦む。操作をミスしてマリオが高所から落下する一瞬、心臓がドキっとする。

まぁマリオにしてみても、こんなヘタクソにプレイされているせいで何度も死ぬ羽目になって、さぞかしいい迷惑だろう。

『はじめてのおつかい』は微笑ましい。

僕は、はじめておつかいをした日の事も、はじめて補助なしの自転車に乗れるようになった日の事も、よく覚えている。

ついでに、はじめて C 言語で printf() 関数を使った日の事も。

相対参照を適当に考慮したHTMLファイルの移動

2012-01-07T19:07:00.000+09:00

1月2日の続きで、自分用ひとこと日記支援システム（？）的なものをだらだら作っています。

ちなみに前回は、「HTML ファイルを生成する仕組みを作ろう」という、なんともややこしい話をしていたのでした。

※ ちなみにあの後、くろねこさんからも非常に有益なアドバイスを頂きました。記事になるところまでできたらここに書きたいと思います。

せっかくですから、ひとこと日記を単に書き溜めるだけでなく、月別のアーカイブページなども生成できれば、情報が整理できて閲覧性も高まるのではないかと考えました（というか、至って平凡な思いつきですが）。

さて、ここで問題になるのが HTML に記述された相対参照の扱いです。

もしも HTML ファイル中に

<img src="baka/aho/img.png" />

とか、

<a href="../../hoge/piyo.html" />

みたいな相対参照が仕組まれていた場合、何も考えずに階層の異なる別ディレクトリ移動・複製すると、リンク切れを起こしてしまうのです。

どうにかしたいですよね。

しかし、この問題を完全に解決しようとすると果てしなく面倒なので、実用上それほど問題が出ないレベルで適当に対処しようというのが本日のテーマです。

【テキストファイルの文字コード判定】

相対リンク云々以前に、テキストファイルを扱う際にまず立ちはだかる第一の関門が文字コードです。テキストに使用される文字コードを適切に判別しないと日本語が化けてしまうのですが、使用されている文字コードを厳密に判断する完全な手段が存在しないのです（たとえば、半角カナのみを特定の配列で使用した場合など）。

ここではとりあえずテキストを読んでみて^※、それを基に文字コードを判定し、それから改めて文字コードを指定してファイルを読みなおす^※※…という無駄な事をやります。

※ テキストの量が多いほど、精度よく文字コードを判定できるため。
※※ 初回の読み込みでは改行文字を無視していましたが、2回目の読み込みでは考慮します。

文字コードを判定するメソッドに関しては、[Javaレビュー]レビューで鍛えるJavaコーディング力その７（文字コードチェック）に有用なサンプルがありましたので、そちらを参考にしました。

String getCharSet(StringBuffer text) 
    throws UnsupportedEncodingException {
  String[] charCodeList = {"SJIS", "EUC-JP", "UTF-8"};
  String str = text.toString();
  for(String code : charCodeList) {
    byte[] bytes = str.getBytes(code);
    if(str.equals(new String(bytes, code))) return code;
  }
  throw new UnsupportedEncodingException();
}

今回は割とメジャーな 3 つの文字コード（Shift JIS、EUC-JP、UTF-8）をサポートします。

上記の getCharSet() メソッドを用いて、文字コードを考慮した HTML ファイルを読み込む処理を書くと、以下のようになります（コードでは省略していますが、IOException が投げられる可能性があります）。

// HTMLファイルを取り敢えず読み込み
File file = new File("hoge.html");
if(!file.exists()) throw new FileNotFoundException();

BufferedReader br   = new BufferedReader(new FileReader(file));
String         line = br.readLine();
StringBuffer   src  = new StringBuffer("");

while(line != null) {
  src.append(line);
  line = br.readLine();
}

// 文字コードを取得
String charset;
try {
  charset = getCharSet(src);
} catch(UnsupportedEncodingException ex) { 
  charset = "UTF-8"; 
}

// 文字コードを指定して再度読み込み
br = new BufferedReader(
       new InputStreamReader(
         new FileInputStream(file), charset));

line = br.readLine();
src.delete(0, src.length());
while(line != null) {
  src.append(line + "\n");
  line = br.readLine();
}
System.out.println(src);

これを走らせると、StringBuffer オブジェクト src に HTML ファイルが読み込まれます。

【HTML から相対参照の抜き出し】

次に、先ほど読み込んだ HTML から相対パスを抽出する処理を考えます。

これも厳密に考えようとするとどうにも面倒なので、正規表現を用いて簡易的に抜く事にしましょう。

とりあえず、img タグなどの src 属性、a タグの href 属性、そして object タグの data 属性に対応するために、以下の正規表現を書きました。

((?:src|href|data)[\\s]*=[\\s]*\")([^\"]+)(\")

この表現は、上記の各属性に対して href="../../hoge/piyo.html" のようにマッチします。

ただし、このままでは絶対パスで表記された URI も拾ってしまうため、さらに除外用のパターンが必要です。

^\\/|^(?:https?|ftp):\\/\\/

これは、http://www.baka.com/hoge.html という文字列のうち、水色で着色した部分にマッチします。また、 /script/piyo.js のような、サーバ名を省略した絶対参照にもマッチします。

本来は、一つのパターンで的確に相対パスのみを抜く事ができればよいのですが、先ほど紹介した

HTML からパスを抜き出すパターン
抜き出したパスから絶対参照を除外するパターン

―― のいずれもが簡易的なものであり、完全に意図した通りには機能しない事^※を考えると、2 段階に分割しておいた方が後々の修正もラクになる気がします。

※ 不完全というのは聞き捨てならんかも知れませんが、よしんば完全に機能するものを作っても、HTML の仕様変更に追従できなければ意味がありません。修正のしやすさを優先しました。

加えて、HTML 特有の問題として、ページ内リンクの扱いを適切に処理する必要があります。たとえば、foo.html 内に、以下のようなタグが含まれていたとしましょう。

<a href="#20120107" />

この場合は、先ほどの正規表現で抜き出したパス（#20120107）の前方に、当該 HTML ファイル名（foo.html）を挿し込み、パス情報を foo.html#20120107 に修正します。

なお、相対パスを一旦絶対パス表記に直しておくと、後々都合がよくなります。

ここで、もしも foo.html の所在が判っていれば、参照先の絶対パス表記は以下のように得る事ができます。

[foo.html の親ディレクトリの絶対パス] + [区切り文字] + [foo.html 内の相対パス]

つまり、foo.html が

C:\Users\Tercel\Documents\html\foo.html

―― に存在して、そこから ..\img\bg.png を相対参照している場合、bg.png の絶対パスは

C:\Users\Tercel\Documents\html\..\img\bg.png

―― となります。"html\..\" のあたりがちょっと冗長ですが、これは一時的なものですのでこのままにしておいても構いません。

ソースから相対パスを抜き出して、新たな基準ディレクトリからのパスに変換する Java コードは以下のようになります。

// ========================================

// HTMLテキストからURIパスを抽出するための正規表現
final String pathRegex = "((?:src|href|data)[\\s]*=[\\s]*\")([^\"]+)(\")";

// URIが絶対パスかどうかを（割と適当に）調べる正規表現
final String absoluteURIRegex = "^\\/|^(?:https?|ftp):\\/\\/";

final Matcher pathMatcher = Pattern.compile(pathRegex, 
               Pattern.CASE_INSENSITIVE).matcher(src);
final Pattern absoluteURLPattern = Pattern.compile(absoluteURIRegex);

// HTMLファイルの相対パスを抽出
while(pathMatcher.find()) {
            
 //String header = pathMatcher.group(1);
 String path  = pathMatcher.group(2);
 //String footer = pathMatcher.group(3);
 if(!absoluteURLPattern.matcher(path).find()) {
  // もし相対パスだったら出力
  System.out.print(path);  
 }
}

ちなみに上記コードには、pathMatcher.group(x); という意味深なコードがありますが、これは正規表現中の丸括弧 ( ) でグループ化された要素を取得するための構文です。

たとえば先ほどの正規表現 ――

((?:src|href|data)[\\s]*=[\\s]*\")([^\"]+)(\")

―― 全体にマッチした文字列のうち、((?:src|href|data)[\\s]*=[\\s]*\") にマッチする部分文字列がグループ1、([^\"]+) にマッチする部分文字列がグループ2、残りの (\") がグループ3となります。

この性質を使によって、以下のようにパスのみを巧く加工する事ができるのです。

String newStr = matcher.group(1) +
hoge(matcher.group(2)) +
matcher.group(3);

あとは、抽出に成功したパスを適切に変換するのみとなりました。

【相対パスの変換】

最後に、任意の基点から見たファイルの相対パスを得る方法を考えます。

ここでは例として、基準ディレクトリ C:\Users\Tercel\Documents\html\aaa\bbb から見た C:\Users\Tercel\Documents\html\xxx\yyy\zzz\index.html の相対パスを得る手続きを一つひとつ確認していく事にします。

まず、両者を並べます。

C:\Users\Tercel\Documents\html\aaa\bbb
C:\Users\Tercel\Documents\html\xxx\yyy\zzz\index.html

次に、名前区切り文字（Windows 環境では "\" 記号、Unix 環境では "/" 記号）で各ファイルパスを分割します。

C:\Users\Tercel\Documents\html\aaa\bbb
C:\Users\Tercel\Documents\html\xxx\yyy\zzz\index.html

区切られたパスを先頭から比較し、共通部分を除去します。

C:\Users\Tercel\Documents\html\aaa\bbb
C:\Users\Tercel\Documents\html\xxx\yyy\zzz\index.html

この時点で残った基点の階層の深さは 2（aaa と bbb）なので、相対パスには親ディレクトリを参照する ".." を 2 つ並べて、対象の残ったパスを続けます。

../../xxx/yy/zzz/index.html

あとは、各階層を Unix デフォルトの区切り文字 "/" で区切れば完成。

../../xxx/yy/zzz/index.html

これをメソッド化すると、以下のようなコードになります。

// ある場所から見たファイルの相対パスを返すよ
//
// 引数1 basisPath  ... 相対パスの基準
// 引数2 targetPath ... 対象ファイルのパス
String getRelativePath(String basisPath, String targetPath) {
  
  File basisFile  = new File(basisPath).getAbsoluteFile();
  File targetFile = new File(targetPath).getAbsoluteFile();
  
  // Java実行環境が使用するファイルの区切り文字の取得
  String separator = File.separator;
  
  // セパレータの正規表現（Windows環境の場合は"\"記号をエスケープ）
  String separatorRegex = separator.replaceAll("\\\\", "\\\\\\\\");
  
  
  // 相対パス表現の基準となるディレクトリ
  File basisDir = basisFile.isDirectory() ? 
    basisFile : basisFile.getParentFile();
  String basisDirPath = basisDir.getAbsolutePath();
  
  // ターゲットとなるディレクトリ
  File targetDir;
  String targetFileName;
  if(targetFile.isDirectory()) {
    targetDir      = targetFile;
    targetFileName = "";
  } else {
    targetDir = targetFile.getParentFile();
    targetFileName = targetFile.getName();
  }
  String targetDirPath  = targetDir.getAbsolutePath();

  // ディレクトリ階層を格納
  String[] basis  = basisDirPath.split(separatorRegex);
  String[] target = targetDirPath.split(separatorRegex);
  
  // 相対パスを生成
  StringBuffer path = new StringBuffer("");
  
  int cIndex;
  int limit = Math.min(target.length, basis.length);
  
  for(cIndex = 0; cIndex < limit; ++cIndex)
    if(!target[cIndex].equals(basis[cIndex])) break;

  for(int i = cIndex; i < basis.length; ++i)
    path.append("../");

  for(int i = cIndex; i < target.length; ++i)
    path.append(target[i] + "/");

  path.append(targetFileName);
  
  return path.toString();
}

このメソッドは、だいたい正しく動作します。

例外的に、

基準となるディレクトリや対象のファイルが存在しない場合
相対パスの基準と対象のファイルがそれぞれ異なるドライブである場合

―― は、期待した結果を得る事ができない場合があります。

また、注意すべき点として、引数に相対パスを与えた場合、（相対参照の）基準はデフォルトで Java 仮想マシンの呼び出し元になります。何を言っているか分からない場合は（ごめんなさい）、両者ともに絶対パスを与えた方が無難でしょう。

問題を知っていて放置するのは気が引けるので、一応この 2 点に関する回避策を述べておきます。

まず前者の場合、基準・あるいは対象ファイルが存在する事を事前に調べる事ができます。なければ一時的にディレクトリを作り、パスを得た後でディレクトリを削除するといった方法で対処できます（ちなみに今回の場合は、基準と対象はどちらも実在のディレクトリであるため、この問題に関してはあまり考える必要はないと思います）。

後者の場合、異なるドライブに対して相対参照する事はできないので、適当な例外を投げて終わるくらいが妥当なのでしょう。それ以上のお節介処理は、「相対パスを得る」という名前から想像できない挙動なので、好ましいとは言えません。

というわけで、HTML ファイルから相対パス抜き、さらにそれを別の基点から見た相対パスに変換するための道具が一応揃いました。

めでたし。

冬休みのてきとうアプリケーションづくり

2012-01-02T22:36:00.002+09:00

【プログラミング初め】

ちょっと前に『作りたい』とか言っていた Web サイトの更新支援ツールの開発をてきとうに始めました。

手始めに前回のネタを GUI 化。

秘密鍵の設定画面

秘密鍵の設定は初回起動時だけ（ただし、変更は可能）。

次回の起動時からは下のような認証画面が表示されるようにしました。 GUI めんどくさい。

秘密鍵の認証画面

FTP アカウントの設定ダイアログボックスも作りました。

FTP アカウントの設定画面
設定した内容は秘密鍵で暗号化されます

で。

ここからちょっと新ネタ（？）です。

ただ独り言をつぶやき続けると、それをソーシャルのかけらもない静的な HTML サイトに投稿してくれるという Twitter を超絶的に劣化させたシステムのメイン画面です。

……うん。

なんだろうねこのインタフェースのデザインセンスの無さは。

とりあえず、ツイート感覚で気軽に Web サイトを更新できれば、放置状態になる事もないかなー……と安直に考えたのです。

しかし、『アプリケーションがある程度完成した時点で満足して結局使わない』というのがいつものパターンなので、たぶん今回もそうなるでしょう。

インタフェースができてくるとテンションが上がるのですが、実は地味なところで壁にぶつかっています。

たとえば、生成される中間ファイル^※の管理方法。

※ これらの中間ファイルを、静的な HTML テンプレートに動的に埋め込んでコンテンツを作ろうという魂胆。

今のところ、ディレクトリを「年」「月」「日」に階層化して、そこにファイルを整理しておこうと考えています。

ただし、階層的管理を実装しようとすると、新規フォルダを作る際の処理オーバヘッドや、ファイル中に埋め込まれた相対パスの扱いが厄介な点など、課題も多くなります。

中間ファイルを階層的に管理する様子

作れば作るほど課題が見えてきますが、もしも自分の中で納得できる解決策が見つかったら、適宜ご紹介していこうかな……と考えています。

明けましておめでとうございます

2012-01-01T18:06:00.000+09:00

旧年中は大変お世話になりました。

おかげさまで、旧年中は 47,000PV を超えるアクセスがありました。

本年は、より一層のご愛顧のほどよろしくお願い申し上げます。

2012年1月1日　たーせる

Javaでできるだけ安全にデータを直列化するためのてきとう備忘録

2011-12-30T14:42:00.007+09:00

現在絶賛放置中の Web サイトを今後の発信基地として再開発したいと思い立った事がきっかけで、最近は Web サイトの簡易更新支援ツールの開発を考えている。

簡単な HTML ジェネレータと FTP クライアントを一緒にしたようなものを作りたい。

特に、FTP クライアントを作る際には、ユーザビリティのためにアカウント情報（ユーザ名とかパスワードとか）をローカルに保存しておき、できるだけワンタッチでサーバにアクセスできる仕様にしたいなと考えている。

データのセーブとロードは至るところで使いそうだし、今日はそこらへんの実装をステップバイステップで書き残しておこうと思う。

あ、この記事のサンプルコードには FTP 接続の実装例は含まれていないので、そういう内容を期待している方にはごめんなさい。あくまでデータのアーカイブに関する話題に絞ってもそれなりの分量になってしまったので……。

【備忘録その1: 直列化】

データのアーカイブに関する最も簡単な方法は、以下のようにしてユーザ名とパスワードをそのまま直列化してしまう事であろう。以下に Java のコードを示す。エラー処理が不十分なのはご容赦。

package com.tercel_tech.sample;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.InputStreamReader;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.io.Serializable;

public class Main {
  private static final String FILENAME = "Account.txt";  // ファイル名

  public static void main(String[] args) {  
    File file = new File(FILENAME);
    
    if(!file.exists()) {
      // 事前に登録したアカウント情報がない場合
      System.out.println("ユーザアカウントが存在しません");
      createUserAccount();
    } else {
      // 事前に登録したアカウント情報がある場合
      System.out.println("ユーザアカウントが存在します");
      loadUserAccount();
    }
  }

  // ==============================
  // ユーザアカウントを作ってファイルに保存するメソッド
  // ==============================
  private static void createUserAccount() {
    System.out.print("ユーザ名を入力して下さい: ");

    // 入力を適当に受け付ける
    String name;
    try {
      name = new BufferedReader(
               new InputStreamReader(
                 System.in)).readLine();
    } catch (IOException ex) {
      name = "名無しさん";
    }

    System.out.print("パスワードを入力して下さい: ");

    // 入力を適当に受け付ける
    String passwd;
    try {
      passwd = new BufferedReader(
                 new InputStreamReader(
                   System.in)).readLine();
    } catch (IOException ex) {
      passwd = "aaaa";
    }

    // ユーザ名とパスワードを基に、アカウント情報を作る
    UserData userData = new UserData(name, passwd);

    // アカウント情報をファイルに保存する
    try {
      new ObjectOutputStream(
        new FileOutputStream(
          FILENAME)).writeObject(userData);
      System.out.println("セーブに成功しました");
    } catch (IOException ex) {
      // 書き込みにミスった
    }
  }

  // ==============================
  // 既に存在するユーザアカウント情報を
  //ファイルからロードするメソッド
  // ==============================
  private static UserData loadUserAccount() {
    try {
      UserData data = (UserData) new ObjectInputStream(
        new FileInputStream(FILENAME)).readObject();

      System.out.println("ロードに成功しました");
      System.out.println("ユーザ名: "   + data.getName());
      System.out.println("パスワード: " + data.getPass());
      return data;
    } catch (Exception ex) {
      return null;
    }
  }
}


// ==============================
// ユーザのログイン情報クラス
// ==============================
class UserData implements Serializable {
  private String name;    // ユーザ名
  private String passwd;  // パスワード

  // コンストラクタ
  // --------------------
  // 引数1: ユーザ名
  // 引数2: パスワード
  public UserData(String name, String passwd) {
    this.name   = name;
    this.passwd = passwd;
  }

  // ゲッタ
  public String getName() { return name;   }
  public String getPass() { return passwd; }
}

まずは、ユーザアカウントが存在しない状態でプログラムを走らせる。

【実行例1】

ユーザアカウントが存在しません
ユーザ名を入力して下さい: tercel
パスワードを入力して下さい: hoge
セーブに成功しました

なお、ユーザ名とパスワードは任意で入力する事ができる。

この後、プログラムを再起動させると、既に存在するアカウント情報を自動的にロードして、その内容を表示してくれる。

【実行例2】

ユーザアカウントが存在します
ロードに成功しました
ユーザ名: tercel
パスワード: hoge

一見するとこれでも充分そうだが、自動生成された Account.txt をメモ帳で開くととんでもない事が分かる。

ユーザ名とパスワードが平文のまま保存されていて、明らかにまずい。常識的に考えて何らかの暗号化を施すべきであろう。

【備忘録その2: 秘密鍵暗号】

そこで、アカウント情報を秘密鍵暗号方式で暗号化する（アルゴリズムは Blowfish）。ちなみに、暗号化に使用した秘密鍵は復号するまで破棄してはいけない。

秘密鍵暗号方式を実装するために、UserData クラスを以下のように書き換える（ただし javax.crypto.Cipher と javax.crypto.spec.SecretKeySpec を import する事）。

この時点で UserData クラスの互換性が無くなっているので、もし Account.txt が残っている場合は削除する。

// ==============================
// ユーザのログイン情報クラス
// ==============================
class UserData implements Serializable {
  private final String SECRET_KEY = "HOGEHOGE_TERCEL";  // 秘密鍵
  private byte[] name;    // ユーザ名
  private byte[] passwd;  // パスワード

  // コンストラクタ
  // --------------------
  // 引数1: ユーザ名
  // 引数2: パスワード
  public UserData(String name, String passwd) {
    try {
      this.name   = encrypt(SECRET_KEY, name);
      this.passwd = encrypt(SECRET_KEY, passwd);
    } catch (Exception ex) {
      System.out.println("暗号化に失敗しました");
      this.name   = null;
      this.passwd = null;
    }
  }

  // ゲッタ
  public String getName() {
    String ret;
    try {
      ret = decrypt(SECRET_KEY, name);
    } catch (Exception ex) {
      System.out.println("復号化に失敗しました");
      ret = null;
    }
    return ret;
  }
  public String getPass() {
    String ret;
    try {
      ret = decrypt(SECRET_KEY, passwd);
    } catch (Exception ex) {
      System.out.println("復号化に失敗しました");
      ret = null;
    }
    return ret;
  }

  // ==============================
  // 暗号化
  // ==============================
  private byte[] encrypt(String key, String text) throws Exception {

    SecretKeySpec sksSpec = new SecretKeySpec(key.getBytes(), "Blowfish");

    Cipher cipher = Cipher.getInstance("Blowfish");

    cipher.init(Cipher.ENCRYPT_MODE, sksSpec);
    byte[] encrypted = cipher.doFinal(text.getBytes());

    return encrypted;
  }

  // ==============================
  // 復号化
  // ==============================
  private String decrypt(String key, byte[] encrypted) throws Exception {
    
    SecretKeySpec sksSpec = new SecretKeySpec(key.getBytes(), "Blowfish");
    Cipher cipher = Cipher.getInstance("Blowfish");

    cipher.init(Cipher.DECRYPT_MODE, sksSpec);
    byte[] decrypted = cipher.doFinal(encrypted);

    return new String(decrypted);
  }
}

実行結果は先程と同じだが、保存されるファイルの内容が異なる。アカウント名とパスワードが傍目には判読できなくなっているのだ。

だが、今度は秘密鍵が平文で保存されている。これでは結局、家の鍵を玄関先に放置しておくようなものだ。

これに関しては、秘密鍵をクラスの外で持てばとりあえずはよさそうだが、どのみちプログラム中で秘密鍵を埋め込む事になるので気持ちが悪い。逆コンパイルされたらイチコロだ。

それよりも、アプリケーションを起動した瞬間にログイン情報の復号化を試みる現在の実装では、無関係な第三者がサーバに繋げてしまう事の方が遙かに問題のような気がする。

【備忘録その3: メッセージダイジェスト】

上記の問題に対処するため、さらに次の方針を採る事にする。

まず、事前にユーザが入力した秘密鍵は、適当なハッシュ関数（ここでは MD5）でメッセージダイジェスト化したハッシュ値として保存しておく。ハッシュ値から元の文字列を復元する事は不可能とされているため、たとえ漏れてもそれほど心配する必要はない。

アカウント情報を復元する際には、ユーザに対して改めて秘密鍵の入力を要求する。鍵を知らない第三者を弾くためである。

システムは入力された鍵をメッセージダイジェスト化し、保存してある鍵のハッシュ値との突き合わせを行う事で同一性のチェックを行う。もし両者が等しければ、入力された鍵の平文を用いてアカウント情報を復号化できる。

この方法だと、プログラムを起動するたびに秘密鍵の入力（認証）を要求する事になるので、少々操作性が悪くなる。しかし、手間の割に守れる情報量は多いので、それなりに割に合っているとは思う。

とりあえず、これを実装した結果を以下に示す。前回の Account.txt が残っている場合は削除する。

まずは UserData クラスから。特に hash メソッドがチャームポイントである。

// ==============================
// ユーザのログイン情報クラス
// ==============================
class UserData implements Serializable {
  //private final String SECRET_KEY = "HOGEHOGE_TERCEL";  // 秘密鍵
  private byte[] secretKey; // 秘密鍵（のハッシュ値）
  private byte[] name;      // ユーザ名
  private byte[] passwd;    // パスワード

  // コンストラクタ
  // --------------------
  // 引数1: 秘密鍵（マスターパスワード）
  // 引数1: ユーザ名
  // 引数2: ログインパスワード
  public UserData(String secretKey, String name, String passwd) {
    try {
      this.secretKey = hash(secretKey);
      this.name      = encrypt(secretKey, name);
      this.passwd    = encrypt(secretKey, passwd);
    } catch (Exception ex) {
      System.err.println("暗号化に失敗しました");
      this.secretKey = null;
      this.name      = null;
      this.passwd    = null;
      System.exit(1);
    }
  }

  // ==============================================
  // 入力された秘密鍵が保存されている秘密鍵と等しいかを
  // ハッシュで比較するよ！
  // ==============================================
  public boolean auth(String secretKey) {
    try {
      byte[] md = hash(secretKey);
      if(md.length != this.secretKey.length) return false;
      for(int i = 0; i < md.length; ++ i) {
        if(md[i] != this.secretKey[i]) return false;
      }
      return true;
    } catch (Exception ex) {
      return false;
    }
  }

  // ===========================================
  // 秘密鍵の平文を渡して名前とパスワードを取得する
  // ===========================================
  public String getName(String key) {
    // 最初に認証
    if(!auth(key)) {
      System.err.println("秘密鍵が正しくありません");
      System.err.println("復号化を行いませんでした");
      System.exit(1);
    }

    String ret;
    try {
      ret = decrypt(key, name);
    } catch (Exception ex) {
      System.err.println("復号化に失敗しました");
      ret = null;
    }
    return ret;
  }

  public String getPass(String key) {
    // 最初に認証
    if(!auth(key)) {
      System.err.println("秘密鍵が正しくありません");
      System.err.println("復号化を行いませんでした");
      System.exit(1);
    }

    String ret;
    try {
      ret = decrypt(key, passwd);
    } catch (Exception ex) {
      System.err.println("復号化に失敗しました");
      ret = null;
    }
    return ret;
  }

  // ==============================
  // メッセージダイジェスト化
  // ==============================
  private byte[] hash(String text) throws Exception {
    MessageDigest md = MessageDigest.getInstance("MD5");
    md.update(text.getBytes());
    return md.digest();
  }

  // ==============================
  // 暗号化
  // ==============================
  private byte[] encrypt(String key, String text) throws Exception {

    SecretKeySpec sksSpec = new SecretKeySpec(key.getBytes(), "Blowfish");

    Cipher cipher = Cipher.getInstance("Blowfish");

    cipher.init(Cipher.ENCRYPT_MODE, sksSpec);
    byte[] encrypted = cipher.doFinal(text.getBytes());

    return encrypted;
  }

  // ==============================
  // 復号化
  // ==============================
  private String decrypt(String key, byte[] encrypted) throws Exception {
    
    SecretKeySpec sksSpec = new SecretKeySpec(key.getBytes(), "Blowfish");
    Cipher cipher = Cipher.getInstance("Blowfish");

    cipher.init(Cipher.DECRYPT_MODE, sksSpec);
    byte[] decrypted = cipher.doFinal(encrypted);

    return new String(decrypted);
  }
}

また、これを使用する Main クラスも以下のように書き換える。

public class Main {
  private static final String FILENAME = "Account.txt";  // ファイル名

  private static String secretKey;

  public static void main(String[] args) {  
    File file = new File(FILENAME);
    
    // 秘密鍵を入力させる
    System.out.print("秘密鍵を入力して下さい: ");
    try {
      secretKey = new BufferedReader(
                    new InputStreamReader(
                      System.in)).readLine();
    } catch (IOException ex) {
      System.err.println("秘密鍵の取得に失敗しました");
      System.exit(1);
    }

    if(!file.exists()) {
      // 事前に登録したアカウント情報がない場合
      System.out.println("ユーザアカウントが存在しません");
      createUserAccount();
    } else {
      // 事前に登録したアカウント情報がある場合
      System.out.println("ユーザアカウントが存在します");
      loadUserAccount();
    }
  }

  // ==============================
  // ユーザアカウントを作ってファイルに保存するメソッド
  // ==============================
  private static void createUserAccount() {
    System.out.print("ユーザ名を入力して下さい: ");

    // 入力を適当に受け付ける
    String name;
    try {
      name = new BufferedReader(
               new InputStreamReader(
                 System.in)).readLine();
    } catch (IOException ex) {
      name = "名無しさん";
    }

    System.out.print("パスワードを入力して下さい: ");

    // 入力を適当に受け付ける
    String passwd;
    try {
      passwd = new BufferedReader(
                 new InputStreamReader(
                   System.in)).readLine();
    } catch (IOException ex) {
      passwd = "aaaa";
    }

    // ユーザ名とパスワードを基に、アカウント情報を作る
    UserData userData = new UserData(secretKey, name, passwd);

    // アカウント情報をファイルに保存する
    try {
      new ObjectOutputStream(
        new FileOutputStream(
          FILENAME)).writeObject(userData);
      System.out.println("セーブに成功しました");
    } catch (IOException ex) {
      // 書き込みにミスった
      System.err.println("セーブに失敗しました");
    }
  }

  // ==============================
  // 既に存在するユーザアカウント情報を
  //ファイルからロードするメソッド
  // ==============================
  private static UserData loadUserAccount() {
    try {
      UserData data = (UserData) new ObjectInputStream(
        new FileInputStream(FILENAME)).readObject();

      System.out.println("ロードに成功しました");
      System.out.println("ユーザ名: "   + data.getName(secretKey));
      System.out.println("パスワード: " + data.getPass(secretKey));
      return data;
    } catch (Exception ex) {
      return null;
    }
  }
}

というか、面倒な人（というか僕）のためにだだ書きしたソースも丸ごと置いておこう。コピペすれば動くと思う。

package com.tercel_tech.sample;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.InputStreamReader;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.io.Serializable;
import java.security.MessageDigest;
import javax.crypto.Cipher;
import javax.crypto.spec.SecretKeySpec;

public class Main {
  private static final String FILENAME = "Account.txt";  // ファイル名

  private static String secretKey;

  public static void main(String[] args) {  
    File file = new File(FILENAME);
    
    // 秘密鍵を入力させる
    System.out.print("秘密鍵を入力して下さい: ");
    try {
      secretKey = new BufferedReader(
                    new InputStreamReader(
                      System.in)).readLine();
    } catch (IOException ex) {
      System.err.println("秘密鍵の取得に失敗しました");
      System.exit(1);
    }

    if(!file.exists()) {
      // 事前に登録したアカウント情報がない場合
      System.out.println("ユーザアカウントが存在しません");
      createUserAccount();
    } else {
      // 事前に登録したアカウント情報がある場合
      System.out.println("ユーザアカウントが存在します");
      loadUserAccount();
    }
  }

  // ==============================
  // ユーザアカウントを作ってファイルに保存するメソッド
  // ==============================
  private static void createUserAccount() {
    System.out.print("ユーザ名を入力して下さい: ");

    // 入力を適当に受け付ける
    String name;
    try {
      name = new BufferedReader(
               new InputStreamReader(
                 System.in)).readLine();
    } catch (IOException ex) {
      name = "名無しさん";
    }

    System.out.print("パスワードを入力して下さい: ");

    // 入力を適当に受け付ける
    String passwd;
    try {
      passwd = new BufferedReader(
                 new InputStreamReader(
                   System.in)).readLine();
    } catch (IOException ex) {
      passwd = "aaaa";
    }

    // ユーザ名とパスワードを基に、アカウント情報を作る
    UserData userData = new UserData(secretKey, name, passwd);

    // アカウント情報をファイルに保存する
    try {
      new ObjectOutputStream(
        new FileOutputStream(
          FILENAME)).writeObject(userData);
      System.out.println("セーブに成功しました");
    } catch (IOException ex) {
      // 書き込みにミスった
      System.err.println("セーブに失敗しました");
    }
  }

  // ==============================
  // 既に存在するユーザアカウント情報を
  //ファイルからロードするメソッド
  // ==============================
  private static UserData loadUserAccount() {
    try {
      UserData data = (UserData) new ObjectInputStream(
        new FileInputStream(FILENAME)).readObject();

      System.out.println("ロードに成功しました");
      System.out.println("ユーザ名: "   + data.getName(secretKey));
      System.out.println("パスワード: " + data.getPass(secretKey));
      return data;
    } catch (Exception ex) {
      return null;
    }
  }
}


// ==============================
// ユーザのログイン情報クラス
// ==============================
class UserData implements Serializable {
  //private final String SECRET_KEY = "HOGEHOGE_TERCEL";  // 秘密鍵
  private byte[] secretKey; // 秘密鍵（のハッシュ値）
  private byte[] name;      // ユーザ名
  private byte[] passwd;    // パスワード

  // コンストラクタ
  // --------------------
  // 引数1: 秘密鍵（マスターパスワード）
  // 引数1: ユーザ名
  // 引数2: ログインパスワード
  public UserData(String secretKey, String name, String passwd) {
    try {
      this.secretKey = hash(secretKey);
      this.name      = encrypt(secretKey, name);
      this.passwd    = encrypt(secretKey, passwd);
    } catch (Exception ex) {
      System.err.println("暗号化に失敗しました");
      this.secretKey = null;
      this.name      = null;
      this.passwd    = null;
      System.exit(1);
    }
  }

  // ==============================================
  // 入力された秘密鍵が保存されている秘密鍵と等しいかを
  // ハッシュで比較するよ！
  // ==============================================
  public boolean auth(String secretKey) {
    try {
      byte[] md = hash(secretKey);
      if(md.length != this.secretKey.length) return false;
      for(int i = 0; i < md.length; ++ i) {
        if(md[i] != this.secretKey[i]) return false;
      }
      return true;
    } catch (Exception ex) {
      return false;
    }
  }

  // ===========================================
  // 秘密鍵の平文を渡して名前とパスワードを取得する
  // ===========================================
  public String getName(String key) {
    // 最初に認証
    if(!auth(key)) {
      System.err.println("秘密鍵が正しくありません");
      System.err.println("復号化を行いませんでした");
      System.exit(1);
    }

    String ret;
    try {
      ret = decrypt(key, name);
    } catch (Exception ex) {
      System.err.println("復号化に失敗しました");
      ret = null;
    }
    return ret;
  }

  public String getPass(String key) {
    // 最初に認証
    if(!auth(key)) {
      System.err.println("秘密鍵が正しくありません");
      System.err.println("復号化を行いませんでした");
      System.exit(1);
    }

    String ret;
    try {
      ret = decrypt(key, passwd);
    } catch (Exception ex) {
      System.err.println("復号化に失敗しました");
      ret = null;
    }
    return ret;
  }

  // ==============================
  // メッセージダイジェスト化
  // ==============================
  private byte[] hash(String text) throws Exception {
    MessageDigest md = MessageDigest.getInstance("MD5");
    md.update(text.getBytes());
    return md.digest();
  }

  // ==============================
  // 暗号化
  // ==============================
  private byte[] encrypt(String key, String text) throws Exception {

    SecretKeySpec sksSpec = new SecretKeySpec(key.getBytes(), "Blowfish");

    Cipher cipher = Cipher.getInstance("Blowfish");

    cipher.init(Cipher.ENCRYPT_MODE, sksSpec);
    byte[] encrypted = cipher.doFinal(text.getBytes());

    return encrypted;
  }

  // ==============================
  // 復号化
  // ==============================
  private String decrypt(String key, byte[] encrypted) throws Exception {
    
    SecretKeySpec sksSpec = new SecretKeySpec(key.getBytes(), "Blowfish");
    Cipher cipher = Cipher.getInstance("Blowfish");

    cipher.init(Cipher.DECRYPT_MODE, sksSpec);
    byte[] decrypted = cipher.doFinal(encrypted);

    return new String(decrypted);
  }
}

これを実行してみよう。まずユーザアカウント情報が存在しない状態から。

【実行例1】

秘密鍵を入力して下さい: HOGEHOGE_TERCEL
ユーザアカウントが存在しません
ユーザ名を入力して下さい: tercel
パスワードを入力して下さい: hoge
セーブに成功しました

次に、この状態でアプリケーションを再起動してみる。

【実行例2】

秘密鍵を入力して下さい: HOGEHOGE_TERCEL
ユーザアカウントが存在します
ロードに成功しました
ユーザ名: tercel
パスワード: hoge

秘密鍵が正しいと、ちゃんと復号に成功する。

ふたたび再起動し、今度は秘密鍵を間違えてみる。

【実行例3】

秘密鍵を入力して下さい: hogehoge_tercel
ユーザアカウントが存在します
ロードに成功しました
秘密鍵が正しくありません
復号化を行いませんでした
Java Result: 1

どのみちロードまではしてしまうわけだが（というかロードしないと秘密鍵の突き合わせが行えないので）、異なる秘密鍵の入力をミスると復号化できないので、セキュリティ的にはそこそこ良さげな気がする。

ちなみに、Account.txt をメモ帳で開くとこんな感じになっている。かろうじて判読できるのはクラス名とフィールドだけだ。

これだけでは解読は困難だろう。

また、プログラムに秘密鍵を直接埋め込んでいるわけでもないので逆コンパイルされてもそれほど痛くもかゆくもない。

というわけで、ひとまずこれでログイン情報を外部に保存する最低限の仕組みが一応できたっぽい。

ふぅ…。

【参考サイト】

近況報告

2011-12-28T19:30:00.003+09:00

先日、研究室の忘年会に参加してきた。

仮配属の3年生に、『興味深い研究テーマとかはあるかな？』と訊いたところ、『たーせるさんのテーマは面白いと思います』と言ってもらえた。社交辞令だとは思うけど、なかなか嬉しいものだ。

現在は、学位論文の執筆・校訂やらプレゼン用のスライドづくりやらで、猫の手も借りたいほど忙しい状態だが、それほど焦る気持ちはない。

同時進行で複数の課題を片付ける事が、昔ほど苦手ではなくなったのはよい事だと思う。

うちの B4 で気がかりなのが、理屈を知らずに人の成果に乗っかる研究姿勢が治らないという事。しかもそれを抜いてしまうと、成果といえるほどの成果が残らないのもまた問題。

既存技術を組み合わせるにしても、誰もが思いつくような組み合わせ方では何の面白みも感じない。『そうきたか！』と思わせる意外性の一つでもあればよいのだが。

先日、ついに炭酸飲料を自作できる夢のような道具を買った。コーラとかが作れる。最初は炭酸を吹きこぼしたり、味の加減を失敗したりと色々あったが、今はなんとか上手にできるようになった。

買い物に行くたびに重たいペットボトルを持って帰らなくてよくなった事と、ペットボトルゴミが減った事は個人的にかなり幸せな事である。

ネットワークスペシャリストに合格したからと言って、その日から何かが変わるわけではない。

今にして思えば、資格自体は別にネットワークである必要はなかった。とにかく“苦手な状態からスタートして、短期間でその分野の最難関資格を取った”という事実が自分の中で大きな自信になった。

今はやりたい事が溢れている状態。

これなら受かる！『ネスペ22 （β版）』

2011-12-21T03:21:00.000+09:00

高度情報処理技術者試験の参考書選びは難しい――。

IT パスポートや基本情報技術者のような万人向けの入門資格では、懇切丁寧な解説やかわいいイラストをふんだんに盛り込んだやさしい参考書が充実しており、受験者は自分にあった教材を選ぶ事ができます。

しかし、試験のレベルが高度化するにつれて、こうしたとっつきやすい参考書が減っていき、よく言えば質実剛健、悪く言えば無味乾燥な参考書ばかりになっていきます。受験層を考えれば致し方のない事かも知れませんが、正直、僕には少々きつい。

というわけで、前置きが長くなりましたが、久しぶりの書評ということで、本日は「本物のネットワークスペシャリストになるための最も詳しい過去問解説と合格のコツ」という副題のついた『ネスペ』シリーズをご紹介したいと思います。

『ネスペ22 β版』

この『ネスペ』シリーズは、まずβ版が（主に通販限定で）発売され、その後、増補改訂を経て正式版が全国の書店に並びます。

数あるネットワークスペシャリストの受験参考書の中でも、『ネスペ』シリーズは試験一回分、それも午後試験のみを徹底的に解説した特殊な構成になっている事が大きな特徴で、表紙には以下の注意書きがあります。

●ご購入時の注意事項
・向き不向きがあります
・ある程度の基礎知識が必要です
・ H22年の午後解説のみです

個人的に読んでみて気になった点は、やはり大手出版社が発行する参考書と比較すると推敲が多少甘く、荒削りな文章が残っている事です。

とはいえ、多くの参考書が前年度版の解説をそのまま使い回しているのに対して、『ネスペ』シリーズはほとんど一から書き直されており、一回の出版にかかる労力がなまなかのものではない事は容易に想像がつきます。

その上、敢えて『β版』を出版するという事は、文章表現の質よりも速報性を重視する姿勢を貫いているのでしょう。

そういう意味では、確かに人を選ぶ本かも知れません。

【「問題」と「設問」の2本立て解説】

『ネスペ』シリーズの特色は、設問とその解答の解説に先立ち、問題そのものの解説がなされる事です。

それはまるで、完成品のジグソーパズルを丁寧に分解しながら各ピースがどう組み合わさっていくのかを一つひとつ確認する作業に似ています。

独学では読解が困難な午後 II 試験の問題解説は特に有用です。読者が負担を感じないように、絶妙なタイミングで解説に移る構成は『ネスペ』シリーズの真骨頂でしょう。

また、問題を読んだ後、設問を理解して正解を導き出す行為は、持っている知識のピースを組み立てる作業に似ています。『ネスペ』シリーズでは、予め問題文をひととおり理解させるため、設問に入っても何を問われているのかが分かりやすい。

つまり、解説を2本立て構成にする事によって、

複雑な問題を読み解く力
知識を組み合わせて解答を作る力

の 2 つを効率的に身に付ける事ができるようになっているように思えます。

集中力が切れて問題文を流し読みしたり、漫然と問題を解いたりする事もなくなるでしょう。

これらは決して派手な受験テクニックではないかも知れませんが、正面突破のための強力な武器にはなります。

【意外と重要、メンタル面のアプローチ】

『ドラゴン桜』という漫画作品の中で、特進コースの担任・桜木が「東大は簡単だ」と生徒に言い聞かせます。それも自信満々に。何度も何度も。

実際、あれを読んで「僕でも東大にいけるかも！」と思ってしまった方もいらっしゃると思います。

同じように、『ネスペ22 β版』には、要所要所に「ネットワークスペシャリスト試験は恐くない」というメッセージが散りばめられています。

暗示にかかりやすい方なら、これを繰り返し読むだけでもネットワークスペシャリストが恐くなくなるから不思議。だんだん「僕でも受かりそう」と思えてきます。

ネットワークスペシャリストは一応難関資格。範囲は広いし難易度は高いしで、いくら勉強しても足りないという考えに囚われがちです。

でも、「午後 II なんて、ただの小問の集まり_{（p.146より）}」という認識があれば、必要以上に懼れる事もなくなります。言われてみれば当たり前だけど言われるまで気付かないような些細な事でも、意識の中に定着すれば非常に大きな力になります。

『ネスペ』シリーズは、こうしたメンタル面のサポートが強く、不安の強い受験者を置き去りにしない配慮が感じられます。試験慣れしている方にとっては少しうるさく感じられるかも知れませんが、これも向き不向きの問題でしょう。

というわけで。

ご紹介した『ネスペ21』『ネスペ22 β版』は、ネットワークスペシャリスト受験の際に最も役に立った2冊でした。

ちなみに『ネスペ21』を初めて手に取った当時（2011年8月）は、正直「何が書いてあるのかさっぱり解らない」という絶望的なレベルでした。勉強を始めてみるとあまりにも役立ちすぎて、ライバルには教えたくないと思うようになったのはここだけの話です。

もし、ネットワークの午後対策で何から手をつけてよいのか判らず、途方に暮れている方がいらっしゃいましたら、書店で（もしくは通販で）チェックしてみてはいかがでしょうか。

ネットワークスペシャリスト試験に合格しました

2011-12-16T12:11:00.002+09:00

応援してくださった皆さん、本当にありがとうございました。

【1月5日追記】

合格証書が届きました。応用情報技術者の証書と併せて。

Processingでかんたん3Dアニメーション

2011-12-03T00:00:00.876+09:00

P5 Advent Calendar 2011 3日目を担当しました。

本日のテーマは【かんたん！ 3Dアニメーション】です。みんな大好き Processing で、こんなものを作ってみました。

段ボール製ロボット“ダンボー”が、ただただ歩き続けるだけのシンプルな作品です。

【どうやって作ったの？】

この作品は、たくさんの直方体パーツを時間の経過とともに少しずつ動かして『歩いているように』見せています。

とはいえ、単純に一つひとつのパーツを無秩序に動かせばよいというわけではありません。パーツ間の幾何構造をきちんと連携させるためにはそれなりの工夫が必要です。

そこでこの記事では、 3 次元アニメーションをそれっぽく見せるためのささやかなテクニックをご紹介していこうかと思います。

記事の前半で技術解説を適当にやって、後半で全ソースコードリストを公開します（ソースだけ欲しい方は記事の最後まで一気にスクロールして下さい）。

技術解説は、3 次元を扱う上で基礎中の基礎となるお話から始まり、座標変換の話題からデータ構造の考察へと話が繋がります。高校レベルの簡単なベクトルや行列の概念を知っている方を想定していますが、本文中には数式が登場しませんので、理系以外の方々にもご理解いただけるかと思います。モーションアニメに限らず、考えを広げるヒントになるかも知れません。

また、OpenProcessing の当該ページからはプロジェクトを一括してダウンロードして頂けます。

【3次元オブジェクトを動かす、ということ】

突然ですが、何があっても形状が変化しないものを、物理用語で『剛体』と言います。

作品中のダンボーはすべて剛体のパーツでできています。

さてさて。剛体の動かし方は以下の 2 種類に分けられます。

並進移動　（物体のすべての頂点が、同じ方向に同じ量だけ移動すること）
回転移動　（物体のすべての頂点が、同じ軸^※のまわりを同じ角度だけ回転すること）

※ 回転軸は必ずしも x, y, z いずれかに一致している必要はなく、任意の（単位）ベクトルを回転軸として設定する事ができます（下図参照）。

視覚化するとこんな感じです。

並進移動

回転移動

繰り返しになりますが、どんなに複雑な移動でも上記の 2 つの移動を合成する事で実現できます。

本質はたったこれだけなのですが、座標軸をどう設定するかによって移動に関する制御のしやすさが全然違ってきます。

次節では、そこらへんの処理を都合よく行うために 2 つの座標系を導入する事にしましょう。

【世界座標とローカル座標】

ざっくり言ってしまえば、『ローカル座標』とは移動する物体に常にくっついていく座標系であり、『世界座標』とは決して変化しない座標系の事です。

解釈としては、“（僕から見て）前後左右”のように極めて局所的で自己中心的な位置表現が『ローカル座標系』、“地球上の緯度経度”のような大局的で一意性のある位置表現が『世界座標系』といったところでしょうか。

したがって、頭や胴体などそれぞれのパーツが独自のローカル座標を持つのに対して、世界座標は一つの 3 次元空間に一つしかありません（下図は頭部パーツのローカル座標と世界座標の対比。面倒なので他のパーツは省略）。

↓
（ t 秒後）
↓

上図からも判る通り、パーツがどれだけ並進や回転移動を繰り返しても、そのローカル座標系における x 軸は常に剛体の右側を向き続けますし、z 軸は剛体の後方向を向き続けます。したがって、ローカル座標系を使えばパーツ同士の相対的な位置関係を容易に定義する事ができるのです（たとえば、胴体が移動しても、右腕は常に胴体から見て“右側”に在り続ける事ができます）。

一方、世界座標系はパーツがどれだけ移動を繰り返しても何の影響も受けませんから、あらゆるパーツの絶対的な位置を定義する唯一の基準として使用する事ができます。

なお、ローカル座標と世界座標は相互変換する事ができます。その理論的な背景は 9 月 27 の日記で簡単に紹介していますので、興味のある方はご覧下さい。

先に述べた剛体の 2 つの移動のうち、並進に関しては translate() メソッドが用意されています。座標系に注意しさえすればさほど難しい事ではありませんので、特に説明は不要でしょう。

ただし、回転移動に関してはさらなる補足が必要かも知れません。次節でご紹介します。

【回転 - ヘディング・ピッチ・バンク】

3 次元剛体のあらゆる回転移動は、ローカル座標系における 3 軸（x、y、z 軸）まわりの回転に分解することができます。逆に言えば、3 軸まわりの回転の合成で、あらゆる回転移動を表す事ができるという事です。これを、『オイラー角表現』と呼びます。

たとえば、先ほど例示した回転移動 ――

―― は、以下のように分解できます（下図中の水色の破線は回転軸）。

↓

Processing では、rotateX()、rotateY()、rotateZ() というメソッドが基本軸まわりの回転移動をサポートしており、それらの合成で大概の回転移動を実現する事が可能です。

注意点としては、たとえ各軸まわりの回転量が同じでも、合成する順序が変わると得られる結果は異なったものになってしまいます。数学的に解釈するならば、行列の積演算が交換法則を持たない事と一緒です。したがって、剛体の回転移動を一意に表すには各軸まわりの回転量だけでなく、それらを合成する際の順序も定義せねばなりません。

各軸まわりの回転を合成する順序には、『ロール・ピッチ・ヨー』や『ヘディング・ピッチ・バンク』など、代表的なものが幾つかあります。この作品では、回転の合成順序をヘディング → ピッチ → バンク（下図参照）で統一しています。

初期状態

↓

ヘディング（ y 軸まわり）

↓

ピッチ（ x 軸まわり）

↓

バンク（ z 軸まわり）

↓

合成完了

先に述べた通りにローカル座標系が定義された場合、ヘディングは y 軸まわりの回転、ピッチは x 軸まわりの回転、そしてバンクは z 軸まわりの回転になります。

【ダンボーのためのデータ構造】

ここまでは、一つのパーツのみに着目して解説を行なってきました。

しかし、ダンボーは複数のパーツの組み合わせでできています。これらの部品はひとまとめにして扱いたいものですが、私たちがよく知っている配列^※ではダンボーの身体構造をうまく表現できません。

※ そもそも配列は、同類のオブジェクト群に番号（順序）をつけて管理するためのデータ構造です。オブジェクトごとに何らかの階層関係がある場合は、必ずしも適切なデータ構造とは言えません。

この作品では 2 種類のデータ構造を用いています。

1 つは各パーツの物理的な親子関係を表す『木構造』で、もう 1 つは『連想配列』です。

木構造における各ノードにはパーツの幾何形状や親子間の相対的な位置関係、および前述の角変位が格納されており、親ノードの位置が変わると子ノードもその影響を受けます（シーングラフに近い構造です）。

胴体をルートとした木構造の概念図: なんかこわくなった

このような構造を導入する利点は、Processing に備わっている行列スタックを巧妙に用いて再帰的にレンダリングを行える事や、一つひとつのパーツの位置を（それらの親ノードにおける）ローカル座標空間から見た位置で定義でき^※、扱いが容易になる事などが挙げられます。

※ 繰り返しになりますが、こうする事で『右腕は胴体から見て常に右側、左腕は常に左側』という相対記述が可能になります。

残念なことに Processing は木構造をサポートしていないため、自力で実装する必要があります。このような場合の常套的なテクニックは、木構造のノードを、子ノードの参照リストを持つ自己参照クラスとして定義することです（Node クラスのソースコード参照）。これにより、ルートノードから全ての子ノードに対する階層的な参照関係が生じますから、木構造の再現は容易であると言えるでしょう。

一方、連想配列は各パーツを階層関係から切り離して（つまり単体で）扱いやすくするために用意しました。たとえば、歩行時の腕の振りをプログラムしたい場合は、腕のみに注目したいですね。そこで、腕パーツをすぐさま参照できるように、パーツの（一意な）名称とパーツのデータを紐付けて連想配列に格納しておくのです。

連想配列の概念図: 識別用の名称から、ただちに対応するパーツを参照できる

連想配列は Processing で用意されている HashMap を利用すれば簡単です。

なお、描画すべきパーツは総て木構造に格納する必要があるのに対して、連想配列には必ずしもすべてのパーツを格納する必要はありません。また、木構造の“ノード”と連想配列の“要素”は、どちらも同じオブジェクトを参照しています。ですので、連想配列経由でアクセスしたオブジェクトに何らかの変更を加えたからといって、対応する木構造のノードと整合しなくなる心配はありません。

【最後のしあげ：足を地に着けるために】

胴体パーツを木構造のルートに設定してしまったため、他のパーツは胴体を中心に動く事になります。このままアニメーションさせるとまるで宙ぶらりんですので、きちんとダンボーの足を地に着けてやらねばなりません。

この作品では、レンダリング直前にダンボーの頂点を世界座標に変換した後、“最も下にある頂点”の値を調べ、それが 0 になるようにモデル全体を鉛直方向に並進させています。

地面にめり込んだダンボー

↓

適切な高さに並進移動

↓

完成

各パーツのローカル座標から世界座標への変換は 9 月 27 の日記で述べた通りです。ただし、各頂点を鉛直軸（y軸）の要素のみで比較するために Comparator をカスタマイズしました（Node クラスのソースコード 99 行目 - 104 行目参照）。

【技術解説あとがき】

Processing Advent Calendar 2011 3日目はいかがでしたでしょうか。

実は、参加表明をしてからというもの、題材選びには大いに悩む日々が続きました。『銀河鉄道を走らせよう』とか『クリスマスイルミネーションを作ろう』とかのアイディアはことごとくボツになりました。

〆切直前にがんばって見つけ出した題材ですので、楽しんで頂けたら非常に嬉しいです。

さて、プログラム自体は一晩で書き上げたものの、ノウハウを文章化する事が殊のほか難しく、やむなく端折ってしまったトピックがいくつかあります。

特に回転に関してはいくら説明しても足りないくらいです。簡単のために本文中では言及を避けましたが、ある角変位を表現する式には一意性がありません。また、オイラー角表現は原理上、ジンバルロックという弱点を抱えています。

この作品を改造する場合は、念のためこうした技術上の問題がある点にご留意ください^※。

※ ジンバルロック問題の解消には四元数という虚数のバケモノみたいなやつを使います。この話題はどんなに駆け足で説明しても長くなってしまうのでまたの機会に。

それでは皆様、Let's Enjoy Processing! ∩( ・ω・)∩

【おまちかねソースコード】

…といきたいところですが、今回の作品はソースコードをコピペしただけでは動きません。

ダンボーのモデルにマッピングするテクスチャが必要です。以下のテクスチャ画像に、それぞれ【d_tab.jpg】【d_face.jpg】【d_coin.jpg】【d_normal.jpg】と名前をつけて、data フォルダに保存して下さい。

d_tab.jpg

d_face.jpg

d_coin.jpg

d_normal.jpg

これで準備ができました。あとは以下のプログラムをすべて Processing IDE にコピー & ペーストすれば、ダンボーが動き出します。たぶん。

【メイン】

final int BACKGROUND_COLOR = 0xFFFFEEBB; // 背景色
Danbo danbo;                             // ダンボー
float time;                              // 経過時間の制御用変数


// ========================================
// 初期化
// ========================================
void setup() {
  size(400, 300, P3D);  
  noStroke();  
  
  danbo = new Danbo(); 
}


// ========================================
// 描画
// ========================================
void draw() {
  // ----------------------------------------
  // ライトが常に正面から当たるように設定
  // ----------------------------------------
  camera(); lights();


  // ----------------------------------------
  // カウンタの更新
  // ----------------------------------------
  float angle = radians(time += .05);
  if(angle >= TWO_PI) angle = time = 0;
  
  
  // ----------------------------------------
  // カメラの設置
  // ----------------------------------------
  camera(700 * sin(angle), -600, -700 * cos(angle), 
         0,                -300,  0, 
         0,                 1,    0);
  
  background(BACKGROUND_COLOR);
  
  
  // ----------------------------------------
  // 歩行モーションの更新
  // ----------------------------------------
  danbo.update();
  
  
  // ----------------------------------------
  // リフレクション（おまけ）
  // ----------------------------------------
  pushMatrix();
  applyMatrix(1.0,  0.0, 0.0, 0.0,
              0.0, -1.0, 0.0, 0.0,
              0.0,  0.0, 1.0, 0.0,
              0.0,  0.0, 0.0, 1.0);
  danbo.render();


  camera();
  pushStyle();
  fill(0xBB << 24 | 0xFFFFFF & BACKGROUND_COLOR);
  rect(0, 0, width, height);
  popStyle();
  popMatrix();

  
  // ----------------------------------------
  // ダンボーの描画
  // ----------------------------------------
  fill(0xFFFFFFFF);
  danbo.render();
}

【Node クラス】

class Node {
  public PVector pos;   // 親からの相対位置
  public PVector org;   // 回転の原点座標
  public PVector size;  // 幅、高さ、奥行き
  public PVector rot;   // 各軸まわりの回転角
  
  private List<Node>          child;        // 子ノードの参照
  private Map<String, PImage> texMap;       // テクスチャ
  private PVector[]           localCoords;  // 頂点のローカル座標
  private PVector[]           globalCoords; // 頂点のグローバル座標
  
  private final int NUM_VERTICES   = 8;  // 頂点数
  private final int NUM_PARTITIONS = 5;  // 分割数
  
  
  // ========================================
  // コンストラクタ
  // ========================================
  public Node() {
    pos          = new PVector(0, 0, 0);
    org          = new PVector(0, 0, 0);
    size         = new PVector(1, 1, 1);
    rot          = new PVector(0, 0, 0);
    
    child        = new ArrayList<Node>();
    texMap       = new HashMap<String, PImage>();
    
    localCoords  = new PVector[NUM_VERTICES];
    globalCoords = new PVector[NUM_VERTICES];
    
    for(int i = 0; i < NUM_VERTICES; ++i) {
      localCoords [i] = new PVector();
      globalCoords[i] = new PVector();
    }
  }
  
  
  // ========================================
  // 子要素追加
  // ========================================
  public void addChild(Node child) {
    this.child.add(child);
  }
  
  
  // ========================================
  // テクスチャ設定
  // ========================================
  public void applyTexture(String name, PImage tex) {
    this.texMap.put(name, tex);
  }
  
  
  // ========================================
  // 最も地面に近い頂点の y 座標値を得る
  // ========================================
  public float getGlobalMaxYCoord() {
    return getGlobalMaxYCoord(Float.MIN_VALUE);
  }
  private float getGlobalMaxYCoord(float maxY) {
    float ret = maxY;
    
    // 箱のローカル頂点をセット
    localCoords[0].set(-.5*size.x, -.5*size.y, -.5*size.z);
    localCoords[1].set(-.5*size.x, -.5*size.y,  .5*size.z);
    localCoords[2].set( .5*size.x, -.5*size.y,  .5*size.z);
    localCoords[3].set( .5*size.x, -.5*size.y, -.5*size.z);
    localCoords[4].set(-.5*size.x,  .5*size.y, -.5*size.z);
    localCoords[5].set(-.5*size.x,  .5*size.y,  .5*size.z);
    localCoords[6].set( .5*size.x,  .5*size.y,  .5*size.z);
    localCoords[7].set( .5*size.x,  .5*size.y, -.5*size.z);    
    
    // ----------------------------------------
    // ワールド行列の計算
    // ----------------------------------------
    pushMatrix();
    
    // 並進
    translate(pos.x, pos.y, pos.z);
    
    // 回転
    translate(org.x, org.y, org.z);
    rotateY(rot.y); rotateX(rot.x); rotateZ(rot.z);
    translate(-org.x, -org.y, -org.z);
    
    // ワールド行列を計算
    PMatrix3D worldMat = (PMatrix3D)getMatrix();
    worldMat.preApply(iViewMat);
    
    // ----------------------------------------
    // ローカル座標から世界座標に変換
    // ----------------------------------------
    for(int i = 0; i < NUM_VERTICES; i++)
      worldMat.mult(localCoords[i], globalCoords[i]);
    
    // ----------------------------------------
    // y座標値で世界座標配列をソート（降順）
    // ----------------------------------------
    Arrays.sort(globalCoords, new java.util.Comparator() {
      public int compare(Object p1, Object p2) {
        return ((PVector)p1).y == ((PVector)p2).y ?  0 :
               ((PVector)p1).y >  ((PVector)p2).y ? -1 : 1;
      }
    });
    // 配列の先頭には最もy座標の大きなベクトルが格納されている
    float tmp = globalCoords[0].y;
    
    if(ret < tmp) ret = tmp;
    
    // ----------------------------------------
    // 再帰的に木をトラバース
    // ----------------------------------------
    for(Node n : child) {
      tmp = n.getGlobalMaxYCoord(ret);
      if (ret < tmp) ret = tmp;
    }
    
    popMatrix();
    return ret;
  }
  
  
  // ========================================
  // ノードのレンダリング
  // ========================================
  void render() {
    pushMatrix();
    
    // 並進
    translate(pos.x, pos.y, pos.z);
    
    // 回転
    translate(org.x, org.y, org.z);
    rotateY(rot.y); rotateX(rot.x); rotateZ(rot.z);
    translate(-org.x, -org.y, -org.z);
    
    float dx = size.x / NUM_PARTITIONS;
    float dz = size.z / NUM_PARTITIONS;
    
    
    // ----------------------------------------
    // 前面の描画
    // ----------------------------------------
    PImage tex = texMap.get("FRONT");
    if(tex == null) tex = texMap.get("DEFAULT");
    
    beginShape(QUAD_STRIP);
    if(tex != null) {
      float texdx = (float)(tex.width-1) / NUM_PARTITIONS;
      texture(tex);
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x - i * dx, -.5 * size.y, -.5 * size.z, texdx * i, 0);
        vertex(.5 * size.x - i * dx,  .5 * size.y, -.5 * size.z, texdx * i, tex.height-1);
      }
    } else {
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x - i * dx, -.5 * size.y, -.5 * size.z);
        vertex(.5 * size.x - i * dx,  .5 * size.y, -.5 * size.z);
      }
    }
    endShape();
    
    
    // ----------------------------------------
    // 背面の描画
    // ----------------------------------------
    tex= texMap.get("BACK");
    if(tex == null) tex = texMap.get("DEFAULT");
    
    beginShape(QUAD_STRIP);
    if(tex != null) {
      float texdx = (float)(tex.width-1) / NUM_PARTITIONS;
      texture(tex);
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(-.5 * size.x + i * dx, -.5 * size.y, .5 * size.z, texdx * i, 0);
        vertex(-.5 * size.x + i * dx,  .5 * size.y, .5 * size.z, texdx * i, tex.height-1);
      }
    } else {
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(-.5 * size.x + i * dx, -.5 * size.y, .5 * size.z);
        vertex(-.5 * size.x + i * dx,  .5 * size.y, .5 * size.z);
      }
    }
    endShape();
    
    
    // ----------------------------------------
    // 天面の描画
    // ----------------------------------------
    tex= texMap.get("TOP");
    if(tex == null) tex = texMap.get("DEFAULT");
    
    beginShape(QUAD_STRIP);
    if(tex != null) {
      float texdx = (float)(tex.width-1) / NUM_PARTITIONS;
      texture(tex);
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x - i * dx, -.5 * size.y,  .5 * size.z, texdx * i, 0);
        vertex(.5 * size.x - i * dx, -.5 * size.y, -.5 * size.z, texdx * i, tex.height-1);
      }
    } else {
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x - i * dx, -.5 * size.y,  .5 * size.z);
        vertex(.5 * size.x - i * dx,  .5 * size.y, -.5 * size.z);
      }
    }
    endShape();
    
    
    // ----------------------------------------
    // 底面の描画
    // ----------------------------------------
    tex= texMap.get("BOTTOM");
    if(tex == null) tex = texMap.get("DEFAULT");
    
    beginShape(QUAD_STRIP);
    if(tex != null) {
      float texdx = (float)(tex.width-1) / NUM_PARTITIONS;
      texture(tex);
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x - i * dx, .5 * size.y,  .5 * size.z, texdx * i, 0);
        vertex(.5 * size.x - i * dx, .5 * size.y, -.5 * size.z, texdx * i, tex.height-1);
      }
    } else {
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x - i * dx, .5 * size.y,  .5 * size.z);
        vertex(.5 * size.x - i * dx, .5 * size.y, -.5 * size.z);
      }
    }
    endShape();
    
    
    // ----------------------------------------
    // 左面の描画
    // ----------------------------------------
    tex= texMap.get("LEFT");
    if(tex == null) tex = texMap.get("DEFAULT");
    
    beginShape(QUAD_STRIP);
    if(tex != null) {
      float texdx = (float)(tex.width-1) / NUM_PARTITIONS;
      texture(tex);
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(-.5 * size.x, -.5 * size.y, -.5 * size.z + i * dz, texdx * i, 0);
        vertex(-.5 * size.x,  .5 * size.y, -.5 * size.z + i * dz, texdx * i, tex.height-1);
      }
    } else {
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(-.5 * size.x, -.5 * size.y, -.5 * size.z + i * dz);
        vertex(-.5 * size.x,  .5 * size.y, -.5 * size.z + i * dz);
      }
    }
    endShape();
    
    
    // ----------------------------------------
    // 右面の描画
    // ----------------------------------------
    tex= texMap.get("RIGHT");
    if(tex == null) tex = texMap.get("DEFAULT");
    
    beginShape(QUAD_STRIP);
    if(tex != null) {
      float texdx = (float)(tex.width-1) / NUM_PARTITIONS;
      texture(tex);
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x, -.5 * size.y, .5 * size.z - i * dz, texdx * i, 0);
        vertex(.5 * size.x,  .5 * size.y, .5 * size.z - i * dz, texdx * i, tex.height-1);
      }
    } else {
      for(int i = 0; i <= NUM_PARTITIONS; ++i) {
        vertex(.5 * size.x, -.5 * size.y, .5 * size.z - i * dz);
        vertex(.5 * size.x,  .5 * size.y, .5 * size.z - i * dz);
      }
    }
    endShape();
    
    
    // ----------------------------------------
    // 再帰的に木をトラバース
    // ----------------------------------------
    for(Node n : child) {
      n.render();
    }
    popMatrix();
  }
}

【Danbo クラス】

Map<String, PImage> imgMap;    // 同じ画像を何度もロードしないようにハッシュで管理
PMatrix3D           iViewMat;  // カメラ逆行列


class Danbo {
  private Map<String, Node> bodyMap;   // ボディパーツ用のハッシュ
  private Node              bodyTree;  // ボディパーツ用の木構造（ルートノード）
  private int               time;      // 経過時間の管理用変数
  
  
  // ========================================
  // コンストラクタ
  // ========================================
  public Danbo() {
    createBody();
    time = 0;
  }
  
  
  // ========================================
  // 歩行モーションの更新
  // ----------------------------------------
  // パーツごとに回転量を再設定しているだけだよ
  // ========================================
  public void update() {
    float angle = radians(++time * 4);
    if(angle >= TWO_PI) angle = time = 0;
    
    float sinAngle = sin(angle);
    float cosAngle = cos(angle);

    // ----------------------------------------
    // 胴の振り
    // ----------------------------------------
    final float BODY_AMPLITUDE = radians(1);  // 振り幅の大きさ（胴を振る角度）

    Node body = bodyMap.get("BODY");
    body.rot.set(0, 0, BODY_AMPLITUDE * cosAngle);
    
    // ----------------------------------------
    // 頭の振り
    // ----------------------------------------
    final float HEAD_AMPLITUDE_X = radians(5);  // 振り幅の大きさ（頭を前後に振る角度）
    final float HEAD_AMPLITUDE_Y = radians(2);  // 　　　〃　　　（ 〃 左右　　〃　　）

    Node head = bodyMap.get("HEAD");
    head.rot.set(HEAD_AMPLITUDE_X * abs(sinAngle), HEAD_AMPLITUDE_Y * sinAngle, 0);

    // ----------------------------------------
    // 腕の振り
    // ---------------------------------------
    final float ARM_AMPLITUDE = radians(40);  // 振り幅の大きさ（腕を振る角度）
    
    Node leftArm = bodyMap.get("LEFTARM");
    leftArm.rot.set(ARM_AMPLITUDE * sinAngle, 0, 0);
    
    Node rightArm = bodyMap.get("RIGHTARM");
    rightArm.rot.set(-ARM_AMPLITUDE * sinAngle, 0, 0);
  
    // ----------------------------------------
    // 脚の動き
    // ----------------------------------------    
    final float LEG_AMPLITUDE = radians(30);  // 振り幅の大きさ（脚を振る角度）
    
    Node leftLeg = bodyMap.get("LEFTLEG");
    leftLeg.rot.set(-LEG_AMPLITUDE * sinAngle, 0, 0);
    
    Node rightLeg = bodyMap.get("RIGHTLEG");
    rightLeg.rot.set(LEG_AMPLITUDE * sinAngle, 0, 0);
  }
  
  
  // ========================================
  // ダンボーのレンダリング
  // ========================================
  public void render() {
    // ----------------------------------------
    // カメラ行列を取得し、その逆行列を得る
    // ----------------------------------------
    iViewMat = (PMatrix3D)getMatrix();
    iViewMat.invert();
    
    // 脚が地面に着くよう、y 方向の並進量を求める
    float shiftY = bodyTree.getGlobalMaxYCoord();
    
    pushMatrix();
    translate(0, -shiftY, 0);
    
    bodyTree.render();
    
    popMatrix();
  }
  
  
  // ----------------------------------------
  
  
  // ========================================
  // 画像ファイルを一括ロードして、imgMapで一元管理します
  // これによって、同じ画像を何度もロードしなくてよくなります
  // ========================================
  private void createImageMap() {
    if(imgMap == null) imgMap = new HashMap<String, PImage>();
    
    if(!imgMap.containsKey("NORMAL")) imgMap.put("NORMAL",loadImage("d_normal.jpg"));
    if(!imgMap.containsKey("COIN"))   imgMap.put("COIN",  loadImage("d_coin.jpg"));
    if(!imgMap.containsKey("TAB"))    imgMap.put("TAB",   loadImage("d_tab.jpg"));
    if(!imgMap.containsKey("FACE"))   imgMap.put("FACE",  loadImage("d_face.jpg"));
  }
  
  
  // ========================================
  // ボディパーツ（幾何情報とマテリアル情報）を、
  // ハッシュ表と木構造に格納します。
  // ハッシュ表は身体部位を高速に検索でき、
  // また、木構造は階層的なアフィン変換を容易に実現できます
  // ========================================
  private void createBody() {
    if(bodyMap == null) bodyMap = new HashMap<String, Node>();
    else bodyMap.clear();
    
    createImageMap();
  
    final float MARGIN = 20;
  
    // ----------------------------------------
    // 胴体
    // ----------------------------------------
    final float BODY_WIDTH  = 130;
    final float BODY_HEIGHT = 155;
    final float BODY_DEPTH  = 130;
    
    // 幾何情報の設定
    Node body = new Node();
    body.size.set(BODY_WIDTH, BODY_HEIGHT, BODY_DEPTH);
    
    // テクスチャの設定
    body.applyTexture("DEFAULT", imgMap.get("NORMAL"));
    body.applyTexture("FRONT",   imgMap.get("COIN"));
    
    
    // ----------------------------------------
    // 頭
    // ----------------------------------------
    final float HEAD_WIDTH  = 275;
    final float HEAD_HEIGHT = 180;
    final float HEAD_DEPTH  = 180;
    
    Node head = new Node();
    head.pos.set (0,          -.5 * (BODY_HEIGHT + HEAD_HEIGHT), 0         );
    head.org.set (0,           .5 * HEAD_HEIGHT,                 0         );
    head.size.set(HEAD_WIDTH,  HEAD_HEIGHT,                      HEAD_DEPTH);
    
    head.applyTexture("TOP",     imgMap.get("TAB"));
    head.applyTexture("FRONT",   imgMap.get("FACE"));
    head.applyTexture("DEFAULT", imgMap.get("NORMAL"));
    
    
    // ----------------------------------------
    // 腕
    // ----------------------------------------
    final float ARM_WIDTH  =  45;
    final float ARM_HEIGHT = 180;
    final float ARM_DEPTH  =  50;
    
    Node leftArm = new Node();
    leftArm.pos.set(-.5 * (BODY_WIDTH + ARM_WIDTH + MARGIN),    // x
                    -.5 * (BODY_HEIGHT - ARM_HEIGHT) + MARGIN,  // y
                      0);                                       // z

    leftArm.org.set (0,         -.5 * ARM_HEIGHT, 0        );
    leftArm.size.set(ARM_WIDTH,   ARM_HEIGHT,     ARM_DEPTH);

    leftArm.applyTexture("DEFAULT", imgMap.get("NORMAL"));


    Node rightArm = new Node();
    rightArm.pos.set( .5 * (BODY_WIDTH + ARM_WIDTH + MARGIN),    // x
                     -.5 * (BODY_HEIGHT - ARM_HEIGHT) + MARGIN,  // y
                       0);                                       // z

    rightArm.org.set (0,         -.5 * ARM_HEIGHT, 0        );
    rightArm.size.set(ARM_WIDTH,  ARM_HEIGHT,      ARM_DEPTH);

    rightArm.applyTexture("DEFAULT", imgMap.get("NORMAL"));
    
    
    // ----------------------------------------
    // 脚
    // ----------------------------------------
    final float LEG_HEIGHT = 130;
    final float LEG_WIDTH  =  50;
    final float LEG_DEPTH  = 100;
    
    Node leftLeg = new Node();
    leftLeg.pos.set(-.5 * (LEG_WIDTH + MARGIN),        // x
                     .5 * (BODY_HEIGHT + LEG_HEIGHT),  // y
                     0);                               // z
                     
    leftLeg.org.set (0,         -.5 * LEG_HEIGHT, 0        );
    leftLeg.size.set(LEG_WIDTH,   LEG_HEIGHT,     LEG_DEPTH);
    
    leftLeg.applyTexture("DEFAULT", imgMap.get("NORMAL"));
  
  
    Node rightLeg = new Node();
    rightLeg.pos.set(.5 * (LEG_WIDTH + MARGIN),        // x
                     .5 * (BODY_HEIGHT + LEG_HEIGHT),  // y
                      0);                              // z
                     
    rightLeg.org.set (0,         -.5 * LEG_HEIGHT, 0        );
    rightLeg.size.set(LEG_WIDTH,   LEG_HEIGHT,     LEG_DEPTH);
    
    rightLeg.applyTexture("DEFAULT", imgMap.get("NORMAL"));
    
    
    // ----------------------------------------
    // スカート（？）
    // ----------------------------------------
    final float SKIRT_HEIGHT = .4 * BODY_DEPTH;
    final float SKIRT_ANGLE1 =  radians(30);
    final float SKIRT_ANGLE2 =  radians(20);
    
    Node backSkirt = new Node();
    
    backSkirt.pos.set( 0,                                 // x
                      .5 * (BODY_HEIGHT + SKIRT_HEIGHT),  // y
                      .5 * BODY_DEPTH);                   // z
                      
    backSkirt.org.set (0,            -.5 * SKIRT_HEIGHT, 0);
    backSkirt.rot.set (SKIRT_ANGLE1,  0,                 0);
    backSkirt.size.set(BODY_WIDTH,    SKIRT_HEIGHT,      1);
    
    backSkirt.applyTexture("DEFAULT", imgMap.get("NORMAL"));
    
    
    Node frontSkirt = new Node();

    frontSkirt.pos.set(  0,                                 // x
                        .5 * (BODY_HEIGHT + SKIRT_HEIGHT),  // y
                       -.5 * BODY_DEPTH);                   // z
                       
    frontSkirt.org.set ( 0,            -.5 * SKIRT_HEIGHT, 0);
    frontSkirt.rot.set (-SKIRT_ANGLE1,  0,                 0);
    frontSkirt.size.set( BODY_WIDTH,    SKIRT_HEIGHT,      1);
    
    frontSkirt.applyTexture("DEFAULT", imgMap.get("NORMAL"));
  
  
    Node leftSkirt = new Node();
    
    leftSkirt.pos.set(-.5 * BODY_WIDTH,                    // x
                       .5 * (BODY_HEIGHT + SKIRT_HEIGHT),  // y
                        0);                                // z
    
    leftSkirt.org.set (0, -.5 * SKIRT_HEIGHT, 0           );
    leftSkirt.rot.set (0,   0,                SKIRT_ANGLE2);
    leftSkirt.size.set(1,   SKIRT_HEIGHT,     BODY_DEPTH  );
    
    leftSkirt.applyTexture("DEFAULT", imgMap.get("NORMAL"));
  
  
    Node rightSkirt = new Node();
    
    rightSkirt.pos.set(.5 * BODY_WIDTH,                    // x
                       .5 * (BODY_HEIGHT + SKIRT_HEIGHT),  // y
                        0);                                // z
                       
    rightSkirt.org.set (0, -.5 * SKIRT_HEIGHT,  0           );
    rightSkirt.rot.set (0,   0,                -SKIRT_ANGLE2);
    rightSkirt.size.set(1,   SKIRT_HEIGHT,      BODY_DEPTH  );
    
    rightSkirt.applyTexture("DEFAULT", imgMap.get("NORMAL"));
    
    
    // ----------------------------------------
    // ハッシュ表の構築
    // ----------------------------------------
    bodyMap.put("HEAD",     head);
    bodyMap.put("BODY",     body);
    bodyMap.put("LEFTARM",  leftArm);
    bodyMap.put("RIGHTARM", rightArm);
    bodyMap.put("LEFTLEG",  leftLeg);
    bodyMap.put("RIGHTLEG", rightLeg);
    
    
    // ----------------------------------------
    // 木の構築
    // ----------------------------------------
    //
    // 　　　　頭
    // 　　　　｜
    // 　　　　胴 ← Root
    // 　　／｜　｜＼ 
    // 右腕　右　左　左腕
    // 　　　脚　足
    // ----------------------------------------
    body.addChild(leftLeg);
    body.addChild(rightLeg);
    body.addChild(head);
    body.addChild(leftArm);
    body.addChild(rightArm);
    bodyTree = body;
  
    // おまけ：スカートのようなもの
    // こちらは別にどうでもいいのでハッシュ表には入れない
    body.addChild(backSkirt);
    body.addChild(frontSkirt);
    body.addChild(leftSkirt);
    body.addChild(rightSkirt);
  }
}

【おまけ】

主催者からの紹介文

モーション、エフェクト、3D、パーティクル、計算幾何、とP5で幅広く活動していらっしゃる @tercel_s さんが、Processing Advent Calendar 2011 に参戦！！ 12/3日担当予定 http://t.co/63T4x3C0 #p5advent2011年12月1日 15:15 via webp5info.com
p5info

うわぁなにこのプレシャーΣ(･Д･ﾉ)ﾉ

C++ で 3 次元 Delaunay 分割

2011-11-21T20:54:00.014+09:00

某所にて、『2 次元の Delaunay 分割は簡単だから別にいらない。むしろ 3 次元 Delaunay 分割のソースコードが欲しい』という要望を間接的に頂いたので、みんなが大好きな C++ で作ってみましたよー (´-ω-`)

OpenProcessing に投稿した 3D Delaunay Triangulation よりも、若干ですがソースコードがブラッシュアップされているような気がします。

今回は 3D という事もあり、結果表示用のプログラムには DirectX 9 を使用しました。最新バージョンを使わなかったのは、Windows XP でも実行できるようにという配慮からです。あ、あくまで環境依存のプログラムは GUI 部分のみであり、アルゴリズム（Delaunay3d.h）は標準 C++ のみで実装しているため、他の環境にもそのまま移植できると思います。

使用方法は昨日の 2 次元 Delaunay 分割とほぼ同じで、以下の通りです（Main.cpp のコールバック関数、および render() 関数に実際のサンプルがあります）。

（※ただし、2 次元と 3 次元のモジュールを共存させようとすると、いたるところで名前が衝突するので注意）

頂点セット（std::set<Tercel::Vector> オブジェクト）を宣言します。
先ほどの頂点セットに、適当な頂点オブジェクトを格納します。
三角形セット（std::set<Tercel::Triangle> オブジェクト）を宣言します。
Tercel::Delaunay3d::getDelaunayTriangles() 関数を呼びます。

このとき、getDelaunayTriangles() の第一引数に頂点セットの『参照』を、第2引数には三角形セットの『ポインタ』を渡します。

例によって、三角形（Tercel::Triangle）は、頂点オブジェクトのポインタを保持していますので、頂点セットから要素を削除した場合（あるいはメモリ上の位置が変わるような事が起きた場合）には、三角形の頂点は不正なポインタになってしまいます。

ご注意ください。

【Delaunay3d.h】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// 3次元Delaunay分割
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


#ifndef ___TERCEL_DELAUNAY3D
#define ___TERCEL_DELAUNAY3D

#include <algorithm>
#include <cfloat>
#include <cmath>

#include <map>
#include <set>

namespace Tercel 
{

  // ポインタの参照先を比較するためのファンクタ
  class PtrComp
  {
  public:
    template <class Type> 
    bool operator()(const Type* a, const Type* b) const
    {
      return *a < *b;
    }
  };


  struct Vector {
    double x, y, z;

    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Vector& v) const 
    {
      return x == v.x && y == v.y && z == v.z;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Vector& v) const 
    {
      return x != v.x ? x < v.x :
             y != v.y ? y < v.y :
                        z < v.z;
    }
  };


  struct Circle 
  {
    Vector center;  // 中心座標
    double radius;  // 半径
  };


  class Triangle {
  public:
    const static size_t NUM_VERTEX = 3;
    const Vector* p[NUM_VERTEX];  // 頂点座標
    
    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Triangle& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>(p[i]);
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());

      for(int i = 0; i < NUM_VERTEX; ++i)
        if(!(p1[i] == p2[i])) return false;
      return true;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Triangle& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>(p[i]);
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());
      for(int i = 0; i < NUM_VERTEX - 1; ++i)
        if(!(*p1[i] == *p2[i])) return *p1[i] < *p2[i];
      return *p1[NUM_VERTEX - 1] < *p2[NUM_VERTEX - 1];
    }
  };


  class Tetrahedron {
  public:
    const static size_t NUM_VERTEX = 4;
    const Vector* p[NUM_VERTEX];

    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Tetrahedron& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>(p[i]);
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());
      for(int i = 0; i < NUM_VERTEX; ++i)
        if(!(p1[i] == p2[i])) return false;
      return true;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Tetrahedron& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>( p[i] );
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());

      for(int i = 0; i < NUM_VERTEX - 1; ++i) 
        if(!(*p1[i] == *p2[i])) return *p1[i] < *p2[i];
      return *p1[NUM_VERTEX - 1] < *p2[NUM_VERTEX - 1];
    }

    // ======================================  
    // 他の四面体と共有点を持つか  
    // ====================================== 
    bool hasCommonPoints(const Tetrahedron& t) const
    {
      for(int i = 0; i < NUM_VERTEX; ++i) 
        for(int j = 0; j < NUM_VERTEX; ++j)
          if(*p[i] == *t.p[j]) return true;
      return false;
    }
  };
  
  class Delaunay3d 
  {
    // ======================================  
    // 型定義
    // ====================================== 
    typedef const std::set<Vector>         ConstVertexSet;
    typedef ConstVertexSet::const_iterator ConstVertexIter;

    typedef std::set<Triangle>             TriangleSet;

    typedef std::set<Tetrahedron>          TetraSet;
    typedef TetraSet::iterator             TetraIter;

    typedef std::map<Tetrahedron, bool>    TetraMap;


  private:
    // ======================================  
    // LU分解による三元一次連立方程式の解法
    // ======================================
    static double lu(double a[3][3], int ip[3]) 
    {
      const int n = 3;
      double weight[n];
   
      for(int k = 0; k < n; ++k)
      {
        ip[k] = k;
        double u = 0;
        for(int j = 0; j < n; ++j) 
        {
          double t = fabs(a[k][j]);
          if (t > u) u = t;
        }
        if (u == 0) return 0;
        weight[k] = 1/u;
      }
      double det = 1;
      for(int k = 0; k < n; ++k) 
      {
        double u = -1;
        int m = 0;
        for(int i = k; i < n; ++i) 
        {
          int ii = ip[i];
          double t = fabs(a[ii][k]) * weight[ii];
          if(t>u) 
          { 
            u = t; 
            m = i;
          }
        }
        int ik = ip[m];
        if (m != k) 
        {
          ip[m] = ip[k]; ip[k] = ik;
          det = -det;
        }
        u = a[ik][k]; det *= u;
        if (u == 0) return 0;
        for (int i = k+1; i < n; ++i) 
        {
          int ii = ip[i];
          double t = (a[ii][k] /= u);
          for(int j = k+1; j < n; ++j) 
            a[ii][j] -= t * a[ik][j];
        }
      }
      return det;
    }

    static void solve(double a[3][3], double b[3], int ip[3], double x[3]) 
    {
      const int n = 3;

      for(int i = 0; i < n; ++i) 
      {
        int ii = ip[i];
        double t = b[ii];
        for (int j = 0; j < i; ++j) 
          t -= a[ii][j] * x[j];
        x[i] = t;
      }
      
      for (int i = n-1; i >= 0; --i) 
      {
        double t = x[i];
        int ii = ip[i];
        for(int j = i+1; j < n; ++j) 
          t -= a[ii][j] * x[j];
        x[i] = t / a[ii][i];
      }
    }

    static double gauss(double a[3][3], double b[3], double x[3]) 
    {
      const int n = 3;
      int ip[n];
      double det = lu(a, ip);
   
      if(det != 0) solve(a, b, ip, x);
      return det;
    }

    // ======================================  
    // 与えられた四面体の外接球を求める
    // ======================================
    static void getCircumcircle(const Tetrahedron& tetra, Circle *dst) 
    {
      const Vector* p0 = tetra.p[0];
      const Vector* p1 = tetra.p[1];
      const Vector* p2 = tetra.p[2];
      const Vector* p3 = tetra.p[3];

      double A[3][3] = {
        {p1->x - p0->x, p1->y - p0->y, p1->z - p0->z},
        {p2->x - p0->x, p2->y - p0->y, p2->z - p0->z},
        {p3->x - p0->x, p3->y - p0->y, p3->z - p0->z}
      };

      double b[3] = {
        0.5 * (p1->x*p1->x - p0->x*p0->x + p1->y*p1->y - p0->y*p0->y + p1->z*p1->z - p0->z*p0->z),
        0.5 * (p2->x*p2->x - p0->x*p0->x + p2->y*p2->y - p0->y*p0->y + p2->z*p2->z - p0->z*p0->z),
        0.5 * (p3->x*p3->x - p0->x*p0->x + p3->y*p3->y - p0->y*p0->y + p3->z*p3->z - p0->z*p0->z)
      };

      double x[3] = {0.0, 0.0, 0.0};

      if(gauss(A, b, x) == 0)
      {
        dst->center.x = dst->center.y = dst->center.z = 0;
        dst->radius   = -1;
      } 
      else 
      {
        dst->center.x = x[0];
        dst->center.y = x[1];
        dst->center.z = x[2];

        double dx = x[0] - p0->x;
        double dy = x[1] - p0->y;
        double dz = x[2] - p0->z;
        dst->radius =  sqrt(dx * dx + dy * dy + dz * dz);
      }
    }
    
    // ======================================  
    // 四面体の重複管理
    // ======================================
    static inline void addElementToRedundanciesMap(TetraMap* pRddcMap,  
        const Tetrahedron& t)   
    {  
      TetraMap::iterator it = pRddcMap->find(t);  
      if(it != pRddcMap->end() && it->second)   
      {  
        // 値を (t, true) から (t, false) に置換  
        pRddcMap->erase(it);  
        pRddcMap->insert(TetraMap::value_type(t, false));  
      } 
      else
        pRddcMap->insert(TetraMap::value_type(t, true));  
    }

  public:
    static void getDelaunayTriangles(ConstVertexSet& pVertexSet, 
        TriangleSet* triangleSet) 
    {

      TetraSet tetraSet;

      // ======================================  
      // 巨大な外部四面体を作る  
      // ======================================
      Vector p0, p1, p2, p3;
      Tetrahedron hugeTetrahedron = {&p0, &p1, &p2, &p3};
      //               !! 注意 !!
      // --------------------------------------
      // hugeTetrahedronの要素はローカル変数への
      // ポインタなので、これらの頂点を共有して
      // いるTetrahedronオブジェクトは、この関数
      // を抜ける前に必ず除去しておく必要がある。
      {
        double maxX, maxY, maxZ;
        double minX, minY, minZ;
        maxX = maxY = maxZ = DBL_MIN;
        minX = minY = minZ = DBL_MAX;
        
        // まず、全ての頂点を包含する球を得る

        // 中心座標を求める
        for(ConstVertexIter it = pVertexSet.begin();
          it != pVertexSet.end(); ++it) 
        {
          if(maxX  < it->x) maxX = it->x; 
          if(it->x < minX ) minX = it->x;
          
          if(maxY  < it->y) maxY = it->y; 
          if(it->y < minY ) minY = it->y;
          
          if(maxZ  < it->z) maxZ = it->z; 
          if(it->z < minZ ) minZ = it->z;
        }
        double centerX = 0.5 * (maxX - minX);
        double centerY = 0.5 * (maxY - minY);
        double centerZ = 0.5 * (maxZ - minZ);

        // 半径を求める
        double dx      = centerX - minX;
        double dy      = centerY - minY;
        double dz      = centerZ - minZ;
        double radius  = sqrt(dx * dx + dy * dy + dz * dz) + 0.1;

        // 4つの頂点
        p0.x = centerX;
        p0.y = centerY + 3.0 * radius;
        p0.z = centerZ;

        p1.x = centerX - 2.0 * sqrt(2.0) * radius;
        p1.y = centerY - radius;
        p1.z = centerZ;

        p2.x = centerX + sqrt(2.0) * radius;
        p2.y = centerY - radius;
        p2.z = centerZ + sqrt(6.0) * radius;

        p3.x = centerX + sqrt(2.0) * radius;
        p3.y = centerY - radius;
        p3.z = centerZ - sqrt(6.0) * radius;
      }

      tetraSet.insert(hugeTetrahedron);

      // --------------------------------------  
      // 点を逐次添加し、反復的に四面体分割を行う  
      // --------------------------------------  
      for(ConstVertexIter vIter = pVertexSet.begin(); 
        vIter != pVertexSet.end(); ++vIter) 
      {
        const Vector* p = &*vIter;

        // --------------------------------------  
        // 追加候補の四面体を保持する一時マップ  
        // --------------------------------------
        TetraMap rddcMap;

        // --------------------------------------  
        // 現在の四面体セットから要素を一つずつ取り出して、    
        // 与えられた点が各々の四面体の外接球の中に含まれるかどうか判定    
        // --------------------------------------  
        for(TetraIter tIter = tetraSet.begin();   
          tIter != tetraSet.end();) 
        {
          // 四面体セットから要素を取りだして…  
          Tetrahedron t = *tIter;

          // 外接球を求める
          Circle c;
          getCircumcircle(t, &c);
          double dx  = c.center.x - p->x;
          double dy  = c.center.y - p->y;
          double dz  = c.center.z - p->z;
          double len = sqrt(dx * dx + dy * dy + dz * dz);

          if(len < c.radius)
          {
            Tetrahedron t0 = {p, t.p[0], t.p[1], t.p[2]};
            addElementToRedundanciesMap(&rddcMap, t0);

            Tetrahedron t1 = {p, t.p[0], t.p[2], t.p[3]};
            addElementToRedundanciesMap(&rddcMap, t1);

            Tetrahedron t2 = {p, t.p[0], t.p[3], t.p[1]};
            addElementToRedundanciesMap(&rddcMap, t2);

            Tetrahedron t3 = {p, t.p[1], t.p[2], t.p[3]};
            addElementToRedundanciesMap(&rddcMap, t3);

            tetraSet.erase(tIter++);
          } else ++tIter;
        }

        for(TetraMap::iterator iter = rddcMap.begin();
          iter != rddcMap.end(); ++iter) 
        {
          if(iter->second) tetraSet.insert(iter->first);
        }
      }

      // --------------------------------------  
      // 最後に、外部三角形の頂点を削除  
      // --------------------------------------  
      for(TetraIter tIter = tetraSet.begin();   
        tIter != tetraSet.end(); )   
      {  
        if(hugeTetrahedron.hasCommonPoints(*tIter))  
           tetraSet.erase(tIter++);  
        else ++tIter;  
      }

      // そして、伝説へ…
      for(TetraIter tIter = tetraSet.begin();   
        tIter != tetraSet.end(); ++tIter) 
      {
        Tetrahedron tetra = *tIter;

        Triangle t0 = {tetra.p[0], tetra.p[1], tetra.p[2]};
        triangleSet->insert(t0);

        Triangle t1 = {tetra.p[0], tetra.p[2], tetra.p[3]};
        triangleSet->insert(t1);

        Triangle t2 = {tetra.p[0], tetra.p[3], tetra.p[1]};
        triangleSet->insert(t2);

        Triangle t3 = {tetra.p[1], tetra.p[2], tetra.p[3]};
        triangleSet->insert(t3);
      }
    }
  };
}

#endif

【Main.cpp】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// 3次元Delaunay分割（表示用GUI）
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
#include <windows.h>
#include <d3d9.h>
#include <d3dx9.h>

#include "Delaunay3d.h"

//------ライブラリの設定とか
#pragma comment(lib, "d3dxof.lib")
#pragma comment(lib, "dxguid.lib")
#pragma comment(lib, "d3dx9.lib")
#pragma comment(lib, "d3d9.lib")
#pragma comment(lib, "winmm.lib")

//------ウィンドウを出すために必要な物あれこれ
const LPCTSTR APP_NAME    = TEXT("DXT2_BY_TERCEL");
const LPCTSTR MUTEX_NAME  = TEXT("MUTEX_DXT2_BY_TERCEL");
const LPCTSTR WINDOW_NAME = TEXT("3D Delaunay Triangulation");

const UINT SCREEN_WIDTH  = 800;
const UINT SCREEN_HEIGHT = 600;
const UINT WINDOW_WIDTH  = (SCREEN_WIDTH + 
              GetSystemMetrics(SM_CXFIXEDFRAME) * 2);
const UINT WINDOW_HEIGHT = (SCREEN_HEIGHT + 
              GetSystemMetrics(SM_CYFIXEDFRAME) * 2 + 
              GetSystemMetrics(SM_CYCAPTION));

// ----------------------------------------
// 頂点集合、およびそれらを基に
// Delaunay 分割された三角形を格納しておくSTLコンテナ
// ----------------------------------------
std::set<Tercel::Vector>   vertices;
std::set<Tercel::Triangle> triangles;


//------Direct3Dに必要なものあれこれ
LPDIRECT3D9 g_pD3D = NULL;
LPDIRECT3DDEVICE9 g_pD3DDev = NULL;
D3DPRESENT_PARAMETERS g_D3DPresentParam;


//------関数プロトタイプ
LRESULT CALLBACK WindowProc(HWND, UINT, WPARAM, LPARAM);
HRESULT InitD3D(HWND);
VOID    Render();

//------FVF
struct ColorPoint
{
  D3DXVECTOR3 vCoord; // 座標
  DWORD   color;  // 色
};

//------エントリポイント
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, 
                   LPSTR lpCmdLine, int nCmdShow)
{
  MSG msg;

  // 多重起動防止
  HANDLE hMutex = CreateMutex(NULL, TRUE, MUTEX_NAME);
  if(GetLastError() == ERROR_ALREADY_EXISTS)
  {
    // 多重起動を検知した際には、
    // 既存のウィンドウを最前面に表示してプログラムを終了する
    HWND existingWnd = FindWindow(APP_NAME, NULL);
    if(existingWnd != NULL)
    {
      if(IsIconic(existingWnd))
        ShowWindowAsync(existingWnd, SW_RESTORE);

      SetForegroundWindow(existingWnd);
    }
    return -1;
  }

  // ウィンドウクラス登録 
  WNDCLASSEX wc;

  wc.cbSize = sizeof(WNDCLASSEX);
  wc.style = CS_HREDRAW | CS_VREDRAW;
  wc.lpfnWndProc = WindowProc;
  wc.cbClsExtra = 0;
  wc.cbWndExtra = 0;
  wc.hInstance = hInstance;
  wc.hIcon = (HICON)LoadImage(NULL,
    MAKEINTRESOURCE(IDI_APPLICATION),
    IMAGE_ICON,
    0,
    0,
    LR_DEFAULTSIZE | LR_SHARED);
  wc.hCursor = (HCURSOR)LoadImage(NULL,
    MAKEINTRESOURCE(IDC_ARROW),
    IMAGE_CURSOR,
    0,
    0,
    LR_DEFAULTSIZE | LR_SHARED);
  wc.hbrBackground = NULL;
  wc.lpszMenuName = NULL;
  wc.lpszClassName = APP_NAME;
  wc.hIconSm = (HICON)LoadImage(NULL,
    MAKEINTRESOURCE(IDI_APPLICATION),
    IMAGE_ICON,
    0,
    0,
    LR_DEFAULTSIZE | LR_SHARED);
  if(!RegisterClassEx(&wc)) return -1;

  // --------------------------------------
  // 固定サイズのウィンドウ作成
  // --------------------------------------
  HWND hWnd = CreateWindow(APP_NAME,
    WINDOW_NAME,
    WS_OVERLAPPEDWINDOW^WS_THICKFRAME^WS_MAXIMIZEBOX, // サイズ変更禁止
    CW_USEDEFAULT,
    CW_USEDEFAULT,
    WINDOW_WIDTH,
    WINDOW_HEIGHT,
    NULL,
    NULL,
    hInstance,
    NULL);
  if(!hWnd) return -1;
  ValidateRect(hWnd, NULL);

  // --------------------------------------
  // Direct3Dを初期化
  // --------------------------------------
  if(FAILED(InitD3D(hWnd))) return -1;

  // --------------------------------------
  // ウィンドウ表示
  // --------------------------------------
  ShowWindow(hWnd, nCmdShow);
  UpdateWindow(hWnd);

  while(TRUE)
  {
    if (PeekMessage(&msg, NULL, 0U, 0U, PM_REMOVE))
    {
      if (msg.message == WM_QUIT) break;
      TranslateMessage(&msg);
      DispatchMessage(&msg);
    }
    else
    {
      Render();
    }
  }
  
  if (g_pD3DDev) g_pD3DDev->Release();
  if (g_pD3D)    g_pD3D->Release();
  CloseHandle(hMutex);
  return (int)msg.wParam;
}

//------Direct3Dを初期化します
HRESULT InitD3D(HWND hWnd)
{
  if((g_pD3D = Direct3DCreate9(D3D_SDK_VERSION)) == NULL)
    return E_FAIL;

  D3DDISPLAYMODE d3ddm;
  if(FAILED(g_pD3D->GetAdapterDisplayMode(D3DADAPTER_DEFAULT, &d3ddm)))
    return E_FAIL;

  ZeroMemory(&g_D3DPresentParam, sizeof(g_D3DPresentParam));
  g_D3DPresentParam.BackBufferCount = 1;
  g_D3DPresentParam.Windowed = TRUE;
  g_D3DPresentParam.BackBufferFormat = D3DFMT_UNKNOWN;
  g_D3DPresentParam.SwapEffect = D3DSWAPEFFECT_DISCARD;
  g_D3DPresentParam.EnableAutoDepthStencil = TRUE;
  g_D3DPresentParam.AutoDepthStencilFormat = D3DFMT_D16;

  if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,
    D3DDEVTYPE_HAL,
    hWnd,
    D3DCREATE_HARDWARE_VERTEXPROCESSING,
    &g_D3DPresentParam,
    &g_pD3DDev)))
  {

    if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,
      D3DDEVTYPE_HAL,
      hWnd,
      D3DCREATE_SOFTWARE_VERTEXPROCESSING,
      &g_D3DPresentParam,
      &g_pD3DDev)))
    {
      if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,
        D3DDEVTYPE_REF,
        hWnd,
        D3DCREATE_SOFTWARE_VERTEXPROCESSING,
        &g_D3DPresentParam,
        &g_pD3DDev)))
      {
        return E_FAIL;
      }
    }
  }

  return S_OK;
}

//------コールバック関数
LRESULT CALLBACK WindowProc(HWND hWnd, UINT uMsg, 
              WPARAM wParam, LPARAM lParam)
{
  switch(uMsg)
  {
    case WM_DESTROY:
      PostQuitMessage(0);
      return 0;

    case WM_CREATE:
      // ======================================  
      // 頂点集合にランダムでデータをセット
      // ======================================  
      for(int i = 0; i < 200; ++i) {
        Tercel::Vector v; // ベクトルを宣言
        float r   = 0.6f * SCREEN_WIDTH + 0.1f * (rand() % 10);
        float phi = D3DXToRadian(rand() % 180 - 90);
        float theta = D3DXToRadian(rand() % 360);

        v.x = r*cos(phi)*cos(theta);
        v.y = r*sin(phi);
        v.z = r*cos(phi)*sin(theta);

        vertices.insert(v);
      }

      // ======================================  
      // Delaunay分割を行う関数
      // ======================================  
      
      // 結果は第2引数（std::set<Tercel::Triangle> 型）に格納される。
      // 
      // Tercel::Triangle のメンバには、3つの頂点p1, p2, p3 があって、
      // いずれも Tercel::Vertex オブジェクトへのポインタになっている
      // これらのポインタは、vertices リストの要素を参照している。
      Tercel::Delaunay3d::getDelaunayTriangles(vertices, &triangles);
      //                                       ~~~~~~~~  ~~~~~~~~~~
      //                                       頂点集合  三角形集合
      //                                       （参照） （ポインタ）
      return 0;
  }
  return DefWindowProc(hWnd, uMsg, wParam, lParam);
}

//------描画関数
VOID Render() {
  if(FAILED(g_pD3DDev->Clear(0, NULL,
    D3DCLEAR_TARGET | D3DCLEAR_ZBUFFER,
    D3DXCOLOR(0.0f, 0.0f, 0.0f, 1.0f),
    1.0f, 0))) return;
  // ワールド行列
  D3DXMATRIX mWorld;
  D3DXMatrixIdentity(&mWorld);
  D3DXMatrixRotationY(&mWorld, 
                      D3DXToRadian(timeGetTime()/50.0f));
  g_pD3DDev->SetTransform(D3DTS_WORLD, &mWorld);

  // ビュー行列
  const float PI = 3.1415926535f;
  D3DXVECTOR3 vEyePoint(0.0f, 0.0f, 1500.0f);
  D3DXVECTOR3 vLookAtPt(0.0f, 0.0f, 0.0f);
  D3DXVECTOR3 vUp(0.0f, 1.0f, 0.0f);

  D3DXMATRIX mView;
  D3DXMatrixLookAtLH(&mView, &vEyePoint, &vLookAtPt, &vUp);
  g_pD3DDev->SetTransform(D3DTS_VIEW, &mView);

  float aspect;
  aspect = (float)SCREEN_WIDTH / (float)SCREEN_HEIGHT;
    
  // 射影行列
  D3DXMATRIX mProj;
  D3DXMatrixPerspectiveFovLH(&mProj,
    D3DXToRadian(45.0f),
    aspect,
    1.0f,
    10000.0f);
  g_pD3DDev->SetTransform(D3DTS_PROJECTION, &mProj);

  if(SUCCEEDED(g_pD3DDev->BeginScene()))
  {
    // ======================================  
    // ここから描画処理
    // ======================================
    ColorPoint vPoint[2];

    vPoint[0].color = D3DCOLOR_ARGB(255, 0, 100, 200);
    vPoint[1].color = D3DCOLOR_ARGB(255, 0, 100, 200);

    g_pD3DDev->SetFVF(D3DFVF_XYZ | D3DFVF_DIFFUSE);
    g_pD3DDev->SetRenderState(D3DRS_LIGHTING, FALSE);

    // ======================================  
    // Delaunay分割された三角形の集合を一つ一つ取りだして描画
    // ======================================  
    for(std::set<Tercel::Triangle>::iterator it = triangles.begin();
      it != triangles.end(); ++it) 
    {
      Tercel::Triangle t = *it;  // 三角形取得
      for(int i = 0; i < 3; ++i) {
        float x1 = (float)t.p[i]->x;
        float y1 = (float)t.p[i]->y;
        float z1 = (float)t.p[i]->z;

        float x2 = (float)t.p[(i+1)%3]->x;
        float y2 = (float)t.p[(i+1)%3]->y;
        float z2 = (float)t.p[(i+1)%3]->z;
        vPoint[0].vCoord = D3DXVECTOR3(x1, y1, z1);
        vPoint[1].vCoord = D3DXVECTOR3(x2, y2, z2);
        g_pD3DDev->DrawPrimitiveUP(D3DPT_LINELIST, 1, 
          vPoint, sizeof(ColorPoint));
      }
    }
    g_pD3DDev->EndScene();
  }

  //----- デバイスロスト対策
  if(g_pD3DDev->Present(NULL, NULL, NULL, NULL) == D3DERR_DEVICELOST)
  {
    if(g_pD3DDev->TestCooperativeLevel() == D3DERR_DEVICELOST)
    {
      return;
    }
    g_pD3DDev->Reset(&g_D3DPresentParam);
  }
}

【おまけ】

勢いでリツイートしてしまったがコード付きという神記事だった2011年11月21日 21:10 via TweetDecknekomatu
nekomatu

わわっΣ(･Д･ﾉ)ﾉ恐縮です（照

続： C++ で Delaunay 分割

2011-11-20T19:36:00.002+09:00

昨日、C++ で Delaunay 分割を実装しました。

で。昨日はデータ構造に線形リスト（std::list）のみを用いており、コーディングしやすい代わりにあまり効率的とはいえない実装になっていました。

今回は、リストの代わりにセット（std::set）を用いる事で、内部処理の効率化を図りました。セットの実体は二分探索木（おそらくは赤黒木）であり、『重複する要素を持たない』『要素の探索が対数時間』というなかなかおいしいデータ構造のため、昨日の段階でネックになっていた冗長な重複判定処理がかなり省かれています。

ではなぜ最初からそれを使わなかったかというと、ずばり std::set の使い方をよく知らなかったからです。……恥ずかしい。

前述の通り、セットの内部構造は二分探索木ですので、セットに格納する要素同士の大小関係を判定できなくてはなりません。単純に比較できない『頂点』や『三角形』の大小関係をどう定義するかで悩みましたが、結局辞書式順序で順序づける事にしました（Tercel::Vector や Tercel::Triangle 構造体の operator< が、比較関数として機能します）。

ソースコードは以下。

【Delaunay2d.h】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// Delaunay分割（やや効率向上？版）
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


#ifndef ___TERCEL_DELAUNAY2D
#define ___TERCEL_DELAUNAY2D

#include <cfloat>
#include <cmath>

#include <map>
#include <set>

namespace Tercel 
{
  struct Vector {
    double x;
    double y;

    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Vector& v) const 
    {
      return (x == v.x && y == v.y);
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Vector& v) const 
    {
      return x != v.x ? x < v.x : y < v.y;
    }
  };


  struct Circle 
  {
    Vector center;  // 中心座標
    double radius;  // 半径
  };


  class Triangle {
  public:
    const Vector* p1, * p2, * p3;  // 頂点座標

  private:
    // 最小の頂点を返す   
    // --------------------------------------
    inline const Vector* getMinVertex() const 
    {
      return *p1 < *p2 && *p1 < *p3 ? p1 : *p2 < *p3 ? p2 : p3;
    }

    // まんなかの頂点を返す
    // --------------------------------------
    inline const Vector* getMidVertex() const 
    {
      return *p1 < *p2 && *p2 < *p3 || *p3 < *p2 && *p2 < *p1 ? p2 : 
             *p2 < *p3 && *p3 < *p1 || *p1 < *p3 && *p3 < *p2 ? p3 : p1;
    }

    // 最大の頂点を返す
    // --------------------------------------
    inline const Vector* getMaxVertex() const 
    {
      return *p2 < *p1 && *p3 < *p1 ? p1 : *p3 < *p2 ? p2 : p3;
    }

  public:
    
    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Triangle& t) const 
    {
      return *p1 == *t.p1 && *p2 == *t.p2 && *p3 == *t.p3 ||
             *p1 == *t.p2 && *p2 == *t.p3 && *p3 == *t.p1 ||
             *p1 == *t.p3 && *p2 == *t.p1 && *p3 == *t.p2 ||

             *p1 == *t.p3 && *p2 == *t.p2 && *p3 == *t.p1 ||
             *p1 == *t.p2 && *p2 == *t.p1 && *p3 == *t.p3 ||
             *p1 == *t.p1 && *p2 == *t.p3 && *p3 == *t.p2;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Triangle& t) const 
    {
      return !(*getMinVertex() == *t.getMinVertex()) ? 
                 *getMinVertex() < *t.getMinVertex() :
             !(*getMidVertex() == *t.getMidVertex()) ? 
                 *getMidVertex() < *t.getMidVertex() :
                 *getMaxVertex() < *t.getMaxVertex();
    }

    // ======================================  
    // 他の三角形と共有点を持つか  
    // ====================================== 
    bool hasCommonPoints(const Triangle& t) const
    {
      return *p1 == *t.p1 || *p1 == *t.p2 || *p1 == *t.p3 ||
             *p2 == *t.p1 || *p2 == *t.p2 || *p2 == *t.p3 ||
             *p3 == *t.p1 || *p3 == *t.p2 || *p3 == *t.p3;
    }
  };

  
  class Delaunay2d 
  {
    // ======================================  
    // 型定義
    // ====================================== 
    typedef const std::set<Vector>         ConstVertexSet;
    typedef ConstVertexSet::const_iterator ConstVertexIter;

    typedef std::set<Triangle>             TriangleSet;
    typedef std::set<Triangle>::iterator   TriangleIter;

    typedef std::map<Triangle, bool>       TriangleMap;

  private:
    // ======================================  
    // 一時マップを使って重複判定  
    // hashMap  
    //  - Key : 三角形  
    //  - Value : 重複していないかどうか  
    //        - 重複していない : true  
    //        - 重複している   : false  
    // ======================================
    static inline void addElementToRedundanciesMap(TriangleMap* pRddcMap,
                                                   Triangle& t) 
    {
      TriangleMap::iterator it = pRddcMap->find(t);
      if(it != pRddcMap->end() && it->second) 
      {
        // 値を (t, true) から (t, false) に置換
        pRddcMap->erase(it);
        pRddcMap->insert(TriangleMap::value_type(t, false));
      } 
      else 
      {
        pRddcMap->insert(TriangleMap::value_type(t, true));
      }
    }

  public:
    static void getDelaunayTriangles(ConstVertexSet& pVertexSet, 
                                     TriangleSet* triangleSet) 
    {

      Triangle hugeTriangle;
      {
        // ======================================  
        // 外部三角形を作る  
        // ======================================  
        double maxX, maxY; maxX = maxY = DBL_MIN;
        double minX, minY; minX = minY = DBL_MAX;
        for(ConstVertexIter it = pVertexSet.begin(); 
            it != pVertexSet.end(); ++it) 
        {
          double x = it->x;
          double y = it->y;
          if(maxX < x) maxX = x; if(minX > x) minX = x;
          if(maxY < y) maxY = y; if(minY > y) minY = y;
        }
        
        // すべての点を包含する矩形の外接円
        double centerX  = (maxX - minX) * 0.5;      // 中心x座標
        double centerY  = (maxY - minY) * 0.5;      // 中心y座標

        double dX   = maxX - centerX;
        double dY   = maxY - centerY;
        double radius = sqrt(dX * dX + dY * dY) + 1.0;  // 半径

        Vector* p1 = new Vector;  // メモリ確保（314行目で解放）
        p1->x    = centerX - sqrt(3.0) * radius;
        p1->y    = centerY - radius;

        Vector* p2 = new Vector;  // メモリ確保（315行目で解放）
        p2->x    = centerX + sqrt(3.0) * radius;
        p2->y    = centerY - radius;

        Vector* p3 = new Vector;  // メモリ確保（316行目で解放）
        p3->x    = centerX;
        p3->y    = centerY + 2.0 * radius;

        hugeTriangle.p1 = p1;
        hugeTriangle.p2 = p2;
        hugeTriangle.p3 = p3;
      }

      triangleSet->insert(hugeTriangle);

      // --------------------------------------  
      // 点を逐次添加し、反復的に三角分割を行う  
      // --------------------------------------  
      for(ConstVertexIter vIter = pVertexSet.begin(); 
          vIter != pVertexSet.end(); ++vIter) 
      {
        const Vector* p = &*vIter;

        // --------------------------------------  
        // 追加候補の三角形を保持する一時マップ
        // -------------------------------------- 
        TriangleMap rddcMap;

        // --------------------------------------  
        // 現在の三角形セットから要素を一つずつ取り出して、  
        // 与えられた点が各々の三角形の外接円の中に含まれるかどうか判定  
        // --------------------------------------  
        for(TriangleIter tIter = triangleSet->begin(); 
            tIter != triangleSet->end();) 
        {
          // 三角形セットから三角形を取りだして…
          Triangle t = *tIter;

          // その外接円を求める。  
          Circle   c;
          {
            // 三角形の各頂点座標を (x1, y1), (x2, y2), (x3, y3) とし、  
            // その外接円の中心座標を (x, y) とすると、  
            //     (x - x1) * (x - x1) + (y - y1) * (y - y1)  
            //   = (x - x2) * (x - x2) + (y - y2) * (y - y2)  
            //   = (x - x3) * (x - x3) + (y - y3) * (y - y3)  
            // より、以下の式が成り立つ  
            //  
            // x = { (y3-y1) * (x2*x2-x1*x1+y2*y2-y1*y1)  
            //     + (y1-y2) * (x3*x3-x1*x1+y3*y3-y1*y1) } / m
            //  
            // y = { (x1-x3) * (x2*x2-x1*x1+y2*y2-y1*y1)  
            //     + (x2-x1) * (x3*x3-x1*x1+y3*y3-y1*y1) } / m  
            //  
            // ただし、  
            // m = 2 * {(x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)} 

            double x1 = t.p1->x;  double y1 = t.p1->y;
            double x2 = t.p2->x;  double y2 = t.p2->y;
            double x3 = t.p3->x;  double y3 = t.p3->y;

            double m = 2.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));  
            double x = ((y3-y1)*(x2*x2-x1*x1+y2*y2-y1*y1)
                       +(y1-y2)*(x3*x3-x1*x1+y3*y3-y1*y1))/m;
            double y = ((x1-x3)*(x2*x2-x1*x1+y2*y2-y1*y1)
                       +(x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1))/m;

            c.center.x = x; 
            c.center.y = y;

            // 外接円の半径 r は、半径から三角形の任意の頂点までの距離に等しい 
            double dx   = t.p1->x - x;
            double dy   = t.p1->y - y;
            double radius = sqrt(dx * dx + dy * dy);

            c.radius = radius;
          }

          double dx = c.center.x - p->x;
          double dy = c.center.y - p->y;
          double dist = sqrt(dx * dx + dy * dy);

          // --------------------------------------  
          // 追加された点が外接円内部に存在する場合、  
          // その外接円を持つ三角形をリストから除外し、  
          // 新たに分割し直す  
          // -------------------------------------- 
          if(dist < c.radius) 
          {
            Triangle t1;
            t1.p1 = p; t1.p2 = t.p1; t1.p3 = t.p2;
            addElementToRedundanciesMap(&rddcMap, t1);
 
            Triangle t2;
            t2.p1 = p; t2.p2 = t.p2; t2.p3 = t.p3;
            addElementToRedundanciesMap(&rddcMap, t2);
 
            Triangle t3;
            t3.p1 = p; t3.p2 = t.p3; t3.p3 = t.p1;
            addElementToRedundanciesMap(&rddcMap, t3);

            triangleSet->erase(tIter++);
          } 
          else ++tIter;
        }

        // --------------------------------------  
        // 一時マップのうち、重複のないものを三角形リストに追加   
        // --------------------------------------  
        for(TriangleMap::iterator iter = rddcMap.begin();
          iter != rddcMap.end(); ++iter) 
        {
          if(iter->second) triangleSet->insert(iter->first);
        }
      }

      // --------------------------------------  
      // 最後に、外部三角形の頂点を削除
      // --------------------------------------  
      for(TriangleIter tIter = triangleSet->begin(); 
          tIter != triangleSet->end(); ) 
      {
        if(hugeTriangle.hasCommonPoints(*tIter))
          triangleSet->erase(tIter++);
        else ++tIter;
      }

      // 巨大三角形の頂点を解放
      delete hugeTriangle.p1;
      delete hugeTriangle.p2;
      delete hugeTriangle.p3;
    }
  };
}

#endif

ちなみに使い方（というか Windows の GUI による結果表示用プログラム）はこちら。

【Main.cpp】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// Delaunay分割（やや効率化版）
//           結果表示用のウィンドウ
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
#include <windows.h>

#include "Delaunay2d.h"

// アプリケーションとウィンドウの初期化に必要な定数
// ----------------------------------------
const LPCTSTR APP_NAME   = TEXT("APP_TERCEL_TECH");
const LPCTSTR MUTEX_NAME = TEXT("MUTEX_TERCEL_TECH");

// クライアント領域の幅と高さ
// ----------------------------------------
const UINT WIDTH  = 400;
const UINT HEIGHT = 300;

// ----------------------------------------
// 頂点集合、およびそれらを基に
// Delaunay 分割された三角形を格納しておくSTLコンテナ
// ----------------------------------------
std::set<Tercel::Vector>   vertices;  // ←型に注意！！
std::set<Tercel::Triangle> triangles; // ←型に注意！！


// コールバック関数（メッセージ処理）
// ========================================
LRESULT CALLBACK WindowProc(HWND hWnd, UINT uMsg,
              WPARAM wParam, LPARAM lParam)
{
  switch(uMsg)
  {
  case WM_DESTROY:
    PostQuitMessage(0);
    return 0;

  case WM_CREATE:
    {
      // ======================================  
      // 頂点集合にランダムでデータをセット
      // ======================================  
      for(int i = 0; i < 200; ++i) {
        Tercel::Vector v; // ベクトルを宣言
        v.x   = rand() % WIDTH;
        v.y   = rand() % HEIGHT;
        vertices.insert(v); // ←変更！！
      }

      // ======================================  
      // Delaunay分割を行う関数
      // ======================================  
      
      // 結果は第2引数（std::set<Tercel::Triangle> 型）に格納される。
      // 
      // Tercel::Triangle のメンバには、3つの頂点p1, p2, p3 があって、
      // いずれも Tercel::Vertex オブジェクトへのポインタになっている
      // これらのポインタは、vertices リストの要素を参照している。
      Tercel::Delaunay2d::getDelaunayTriangles(vertices, &triangles);
      //                     ~~~~~~~~  ~~~~~~~~~~
      //                     頂点集合  三角形集合
      //                     （参照） （ポインタ）
    }
    return 0;

  case WM_PAINT:
    {
      typedef std::set<Tercel::Triangle> TriangleSet;

      PAINTSTRUCT ps;
      HDC     hdc  = BeginPaint(hWnd, &ps);
      HPEN    hPen = CreatePen(PS_SOLID, 1, RGB(0, 0, 0));
      
      SelectObject(hdc, hPen);

      // ======================================  
      // Delaunay分割された三角形の集合を一つ一つ描画
      // ======================================  
      for(TriangleSet::iterator it = triangles.begin();
        it != triangles.end(); ++it) 
      {
        Tercel::Triangle t = *it;  // 三角形取得
        int x, y;

        x = int(t.p1->x + 0.5);
        y = int(t.p1->y + 0.5);

        MoveToEx(hdc, x, y, NULL);

        x = int(t.p2->x + 0.5);
        y = int(t.p2->y + 0.5);

        LineTo(hdc, x, y);
        MoveToEx(hdc, x, y, NULL);

        x = int(t.p3->x + 0.5);
        y = int(t.p3->y + 0.5);

        LineTo(hdc, x, y);
        MoveToEx(hdc, x, y, NULL);

        x = int(t.p1->x + 0.5);
        y = int(t.p1->y + 0.5);

        LineTo(hdc, x, y);
      }
      DeleteObject(hPen);
      EndPaint(hWnd, &ps);      
    }
    return 0;
  }
  return DefWindowProc(hWnd, uMsg, wParam, lParam);
}


// エントリポイント
// ========================================
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,
           LPSTR lpCmdLine, int nCmdShow)
{

  // 多重起動防止
  HANDLE hMutex = CreateMutex(NULL, TRUE, MUTEX_NAME);
  if(GetLastError() == ERROR_ALREADY_EXISTS)
  {
    // 多重起動を検知した際には、
    // 既存のウィンドウを最前面に表示してプログラムを終了する
    HWND existingWnd = FindWindow(APP_NAME, NULL);
    if(existingWnd != NULL)
    {
      if(IsIconic(existingWnd))
        ShowWindowAsync(existingWnd, SW_RESTORE);

      SetForegroundWindow(existingWnd);
    }
    return -1;
  }

  // ウィンドウクラス登録
  WNDCLASSEX wc;
  wc.cbSize    = sizeof(WNDCLASSEX);
  wc.style     = CS_HREDRAW | CS_VREDRAW;
  wc.lpfnWndProc   = WindowProc;
  wc.cbClsExtra  = 0;
  wc.cbWndExtra  = sizeof(LONG_PTR);
  wc.hInstance   = hInstance;
  wc.hIcon     = (HICON)LoadImage(NULL,
               MAKEINTRESOURCE(IDI_APPLICATION),
               IMAGE_ICON,
               0,
               0,
               LR_DEFAULTSIZE | LR_SHARED);
  wc.hCursor     = (HICON)LoadImage(NULL,
               MAKEINTRESOURCE(IDC_ARROW),
               IMAGE_CURSOR,
               0,
               0,
               LR_DEFAULTSIZE | LR_SHARED);
  wc.hbrBackground = (HBRUSH)COLOR_BACKGROUND + 1;
  wc.lpszMenuName  = NULL;
  wc.lpszClassName = APP_NAME;
  wc.hIconSm = (HICON)LoadImage(NULL,
               MAKEINTRESOURCE(IDI_APPLICATION),
               IMAGE_ICON,
               0,
               0,
               LR_DEFAULTSIZE | LR_SHARED);
  if(!RegisterClassEx(&wc)) return -1;

  // クライアント領域のサイズから、ウィンドウそのものの幅と高さを計算
  UINT windowWidth  = WIDTH  + GetSystemMetrics(SM_CXFIXEDFRAME) * 2;
  UINT windowHeight = HEIGHT + GetSystemMetrics(SM_CYFIXEDFRAME) * 2
                 + GetSystemMetrics(SM_CYCAPTION);
  // ウィンドウ生成
  HWND hWnd = CreateWindow(APP_NAME,
          TEXT("はじめてのDelaunay分割"),
          WS_OVERLAPPEDWINDOW^WS_THICKFRAME^WS_MAXIMIZEBOX,
          CW_USEDEFAULT,
          CW_USEDEFAULT,
          windowWidth,
          windowHeight,
          NULL,
          NULL,
          hInstance,
          NULL);
  if(!hWnd) return -1;

  // ウィンドウ表示
  ShowWindow(hWnd, nCmdShow);
  UpdateWindow(hWnd);

  // メッセージループ
  MSG msg;
  while(TRUE)
  {
    if(PeekMessage(&msg, NULL, 0U, 0U, PM_REMOVE))
    {
      if(msg.message == WM_QUIT) break;
      TranslateMessage(&msg);
      DispatchMessage(&msg);
    }
  }
  CloseHandle(hMutex);

  return (int)msg.wParam;
}

イテレータのお陰で、データ構造が変わっても Main.cpp 側のソースコードはほとんどそのままです。便利ですね。

C++ で Delaunay 分割（ただし2次元）

2011-11-19T18:11:00.003+09:00

ちょっと学業の方が忙しくなってしまい、ウェブログを放置していました……。

忘れていたわけじゃないんだよ。ただちょっとモチベーションが足りなかっただけ。

今日は、某所（CG 系の研究室）で需要がありそうだったので、以前 Processing で書いていた Delaunay 分割を C++ で実装し直してみる事にしました。

ただし、C++ はものすごく苦手なので、ちょっと残念なソースになってしまっています。 ⇒ 11月20日追記：もうちょっと内部実装がマシになりました。

ソースコード上部の「view plain」という文字をクリックするとコピペ用の窓が出ますので、使いたい方はご自由にどうぞ。

【アルゴリズム部分（Delaunay2d.h）】

#ifndef ___TERCEL_DELAUNAY2D
#define ___TERCEL_DELAUNAY2D

#include <cfloat>
#include <cmath>
#include <list>

namespace Tercel 
{
    struct Vector 
    {
        double x, y;

        bool operator==(const Vector& v) const 
        {
            return (x == v.x && y == v.y);
        }
    };


    struct Circle 
    {
        Vector center;  // 中心座標
        double radius;  // 半径
    };


    class Triangle 
    {
    public:
        const Vector* p1, * p2, * p3;    // 頂点座標

    public:
        // ======================================  
        // 等価性の判定
        // ======================================   
        bool operator==(const Triangle& t) const 
        {
            return(*p1 == *t.p1 && *p2 == *t.p2 && *p3 == *t.p3 ||
                   *p1 == *t.p2 && *p2 == *t.p3 && *p3 == *t.p1 ||
                   *p1 == *t.p3 && *p2 == *t.p1 && *p3 == *t.p2 ||

                   *p1 == *t.p3 && *p2 == *t.p2 && *p3 == *t.p1 ||
                   *p1 == *t.p2 && *p2 == *t.p1 && *p3 == *t.p3 ||
                   *p1 == *t.p1 && *p2 == *t.p3 && *p3 == *t.p2 );
        }

        // ======================================  
        // 他の三角形と共有点を持つか  
        // ======================================   
        bool hasCommonPoints(const Triangle& t) 
        {
            return(*p1 == *t.p1 || *p1 == *t.p2 || *p1 == *t.p3 ||
                   *p2 == *t.p1 || *p2 == *t.p2 || *p2 == *t.p3 ||
                   *p3 == *t.p1 || *p3 == *t.p2 || *p3 == *t.p3 );
        }
    };

    
    class Delaunay2d 
    {
    private:
        static void manageDuplicativeTriangles(std::list<Triangle>* newTriangleList,
                                        std::list<Triangle>* duplicativeTriangleList, 
                                        const Triangle& t) 
        {
            typedef std::list<Triangle>::iterator triIter;

            bool existsInNewTriangleList = false;
            for(triIter iter = newTriangleList->begin();
                iter != newTriangleList->end(); iter++) 
            {
                
                if(*iter == t) {
                    existsInNewTriangleList = true;
                    bool existsInDuplicativeTriangleList = false;

                    for(triIter iter2 = duplicativeTriangleList->begin();
                        iter2 != duplicativeTriangleList->end(); iter2++) 
                    {
                        if(*iter2 == t) 
                        { 
                            existsInDuplicativeTriangleList = true;
                            break;
                        }

                    }
                    if(!existsInDuplicativeTriangleList) 
                    {
                        duplicativeTriangleList->push_back(t);
                    }
                    break;
                }
            }
            if(!existsInNewTriangleList) newTriangleList->push_back(t);
        }

    public:
        static void getDelaunayTriangles(const std::list<Vector>& vertexList, 
            std::list<Triangle>* triangleList) 
        {
            typedef std::list<Vector>::const_iterator cVtxIter;
            typedef std::list<Triangle>::iterator     triIter;

            Triangle hugeTriangle;
            {
                // ======================================  
                // 外部三角形を作る  
                // ======================================  
                double maxX, maxY; maxX = maxY = DBL_MIN;
                double minX, minY; minX = minY = DBL_MAX;
                for(cVtxIter it = vertexList.begin(); it != vertexList.end(); ++it) 
                {
                    double x = it->x;
                    double y = it->y;
                    if(maxX < x) maxX = x; if(minX > x) minX = x;
                    if(maxY < y) maxY = y; if(minY > y) minY = y;
                }
                
                // すべての点を包含する矩形の外接円
                double centerX  = (maxX - minX) * 0.5;          // 中心x座標
                double centerY  = (maxY - minY) * 0.5;          // 中心y座標

                double dX     = maxX - centerX;
                double dY     = maxY - centerY;
                double radius = sqrt(dX * dX + dY * dY) + 1.0;  // 半径

                Vector* p1 = new Vector;    // メモリ確保（266行目で解放）
                p1->x      = centerX - sqrt(3.0) * radius;
                p1->y      = centerY - radius;

                Vector* p2 = new Vector;    // メモリ確保（267行目で解放）
                p2->x      = centerX + sqrt(3.0) * radius;
                p2->y      = centerY - radius;

                Vector* p3 = new Vector;    // メモリ確保（268行目で解放）
                p3->x      = centerX;
                p3->y      = centerY + 2.0 * radius;

                hugeTriangle.p1 = p1;
                hugeTriangle.p2 = p2;
                hugeTriangle.p3 = p3;
            }

            triangleList->push_back(hugeTriangle);

            // --------------------------------------  
            // 点を逐次添加し、反復的に三角分割を行う  
            // --------------------------------------  
            for(cVtxIter vIter = vertexList.begin(); vIter != vertexList.end(); ++vIter) 
            {
                const Vector* p = &(*vIter);

                // --------------------------------------  
                // 追加候補の三角形を保持する一時リスト  
                // --------------------------------------  
                std::list<Triangle> newTriangleList;          // 新規分割された三角形                
                std::list<Triangle> duplicativeTriangleList;  // 重複リスト
                

                // --------------------------------------  
                // 現在の三角形セットから要素を一つずつ取り出して、  
                // 与えられた点が各々の三角形の外接円の中に含まれるかどうか判定  
                // --------------------------------------  
                for(triIter tIter = triangleList->begin(); tIter != triangleList->end(); ) 
                {
                    // 三角形セットから三角形を取りだして…
                    Triangle t = *tIter;

                    // その外接円を求める。  
                    Circle   c;
                    {
                        // 三角形の各頂点座標を (x1, y1), (x2, y2), (x3, y3) とし、  
                        // その外接円の中心座標を (x, y) とすると、  
                        //     (x - x1) * (x - x1) + (y - y1) * (y - y1)  
                        //   = (x - x2) * (x - x2) + (y - y2) * (y - y2)  
                        //   = (x - x3) * (x - x3) + (y - y3) * (y - y3)  
                        // より、以下の式が成り立つ  
                        //  
                        // x = { (y3 - y1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)  
                        //     + (y1 - y2) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)} / c  
                        //  
                        // y = { (x1 - x3) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)  
                        //     + (x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)} / c  
                        //  
                        // ただし、  
                        //   c = 2 * {(x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)} 

                        double x1 = t.p1->x;    double y1 = t.p1->y;
                        double x2 = t.p2->x;    double y2 = t.p2->y;
                        double x3 = t.p3->x;    double y3 = t.p3->y;

                        double m = 2.0 * ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1));  
                        double x = ((y3 - y1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)  
                                  + (y1 - y2) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)) / m;  
                        double y = ((x1 - x3) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)  
                                  + (x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)) / m;

                        c.center.x = x; 
                        c.center.y = y;

                        // 外接円の半径 r は、半径から三角形の任意の頂点までの距離に等しい 
                        double dx     = t.p1->x - x;
                        double dy     = t.p1->y - y;
                        double radius = sqrt(dx * dx + dy * dy);

                        c.radius = radius;
                    }

                    double dx   = c.center.x - p->x;
                    double dy   = c.center.y - p->y;
                    double dist = sqrt(dx * dx + dy * dy);

                    // ======================================  
                    // 一時リストを使って重複判定  
                    // ======================================  
                    if(dist < c.radius) 
                    {
                        // 再分割

                        Triangle t1;
                        t1.p1 = p; t1.p2 = t.p1; t1.p3 = t.p2;
                        manageDuplicativeTriangles(&newTriangleList, &duplicativeTriangleList, t1);

                        Triangle t2;
                        t2.p1 = p; t2.p2 = t.p2; t2.p3 = t.p3;
                        manageDuplicativeTriangles(&newTriangleList, &duplicativeTriangleList, t2);

                        Triangle t3;
                        t3.p1 = p; t3.p2 = t.p3; t3.p3 = t.p1;
                        manageDuplicativeTriangles(&newTriangleList, &duplicativeTriangleList, t3);

                        tIter = triangleList->erase(tIter);
                    } 
                    else ++tIter;
                }

                // --------------------------------------  
                // 一時リストのうち、重複のないものを三角形リストに追加   
                // --------------------------------------  
                for(triIter iter = newTriangleList.begin();
                    iter != newTriangleList.end(); ++iter) 
                {
                    bool exists = false;
                    for(triIter iter2 = duplicativeTriangleList.begin();
                        iter2 != duplicativeTriangleList.end(); ++iter2)
                    {
                        if(*iter == *iter2) 
                        {
                            exists = true;
                            break;
                        }
                    }
                    if(!exists) triangleList->push_back(*iter);
                }
            }

            // 最後に、外部三角形の頂点を削除
            for(triIter tIter = triangleList->begin(); tIter != triangleList->end(); ) 
            {
                if(hugeTriangle.hasCommonPoints(*tIter)) tIter = triangleList->erase(tIter);
                else ++tIter;
            }

            // 巨大三角形の頂点を解放
            delete hugeTriangle.p1;
            delete hugeTriangle.p2;
            delete hugeTriangle.p3;
        }
    };
}

#endif

【テスト用の GUI （Main.cpp）】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// 結果表示用のウィンドウ
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
#include <windows.h>

#include "Delaunay2d.h"

// アプリケーションとウィンドウの初期化に必要な定数
// ----------------------------------------
const LPCTSTR APP_NAME   = TEXT("APP_TERCEL_TECH");
const LPCTSTR MUTEX_NAME = TEXT("MUTEX_TERCEL_TECH");

// クライアント領域の幅と高さ
// ----------------------------------------
const UINT WIDTH  = 400;
const UINT HEIGHT = 300;

// ----------------------------------------
// 頂点集合、およびそれらを基に
// Delaunay 分割された三角形を格納しておくSTLコンテナ
// ----------------------------------------
std::list<Tercel::Vector>   vertices;
std::list<Tercel::Triangle> triangles;


// コールバック関数（メッセージ処理）
// ========================================
LRESULT CALLBACK WindowProc(HWND hWnd, UINT uMsg,
                            WPARAM wParam, LPARAM lParam)
{
    switch(uMsg)
    {
    case WM_DESTROY:
        PostQuitMessage(0);
        return 0;

    case WM_CREATE:
        {
            // ======================================  
            // 頂点集合にランダムでデータをセット
            // ======================================  
            for(int i = 0; i < 200; ++i) {
                Tercel::Vector v;  // ベクトルを宣言
                v.x   = rand() % WIDTH;
                v.y   = rand() % HEIGHT;
                vertices.push_back(v);
            }

            // ======================================  
            // Delaunay分割を行う関数
            // ======================================  
            
            // 結果は第2引数（std::list<Tercel::Triangle> 型）に格納される。
            // 
            // Tercel::Triangle のメンバには、3つの頂点p1, p2, p3 があって、
            // いずれも Tercel::Vertex オブジェクトへのポインタになっている
            // これらのポインタは、vertices リストの要素を参照している。
            Tercel::Delaunay2d::getDelaunayTriangles(vertices, &triangles);
            //                                       ~~~~~~~~  ~~~~~~~~~~
            //                                       頂点集合  三角形集合
            //                                       （参照） （ポインタ）
        }
        return 0;

    case WM_PAINT:
        {
            typedef std::list<Tercel::Triangle> triangleList;

            PAINTSTRUCT ps;
            HDC         hdc  = BeginPaint(hWnd, &ps);
            HPEN        hPen = CreatePen(PS_SOLID, 1, RGB(0, 0, 0));
            
            SelectObject(hdc, hPen);

            // ======================================  
            // Delaunay分割された三角形の集合を一つ一つ描画
            // ======================================  
            for(triangleList::iterator it = triangles.begin();
                it != triangles.end(); ++it) 
            {
                Tercel::Triangle t = *it;  // 三角形取得
                int x, y;

                x = int(t.p1->x + 0.5);
                y = int(t.p1->y + 0.5);

                MoveToEx(hdc, x, y, NULL);

                x = int(t.p2->x + 0.5);
                y = int(t.p2->y + 0.5);

                LineTo(hdc, x, y);
                MoveToEx(hdc, x, y, NULL);

                x = int(t.p3->x + 0.5);
                y = int(t.p3->y + 0.5);

                LineTo(hdc, x, y);
                MoveToEx(hdc, x, y, NULL);

                x = int(t.p1->x + 0.5);
                y = int(t.p1->y + 0.5);

                LineTo(hdc, x, y);
            }
            DeleteObject(hPen);
            EndPaint(hWnd, &ps);            
        }
        return 0;
    }
    return DefWindowProc(hWnd, uMsg, wParam, lParam);
}


// エントリポイント
// ========================================
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,
                   LPSTR lpCmdLine, int nCmdShow)
{

    // 多重起動防止
    HANDLE hMutex = CreateMutex(NULL, TRUE, MUTEX_NAME);
    if(GetLastError() == ERROR_ALREADY_EXISTS)
    {
        // 多重起動を検知した際には、
        // 既存のウィンドウを最前面に表示してプログラムを終了する
        HWND existingWnd = FindWindow(APP_NAME, NULL);
        if(existingWnd != NULL)
        {
            if(IsIconic(existingWnd))
                ShowWindowAsync(existingWnd, SW_RESTORE);

            SetForegroundWindow(existingWnd);
        }
        return -1;
    }

    // ウィンドウクラス登録
    WNDCLASSEX wc;
    wc.cbSize        = sizeof(WNDCLASSEX);
    wc.style         = CS_HREDRAW | CS_VREDRAW;
    wc.lpfnWndProc   = WindowProc;
    wc.cbClsExtra    = 0;
    wc.cbWndExtra    = sizeof(LONG_PTR);
    wc.hInstance     = hInstance;
    wc.hIcon         = (HICON)LoadImage(NULL,
                           MAKEINTRESOURCE(IDI_APPLICATION),
                           IMAGE_ICON,
                           0,
                           0,
                           LR_DEFAULTSIZE | LR_SHARED);
    wc.hCursor       = (HICON)LoadImage(NULL,
                           MAKEINTRESOURCE(IDC_ARROW),
                           IMAGE_CURSOR,
                           0,
                           0,
                           LR_DEFAULTSIZE | LR_SHARED);
    wc.hbrBackground = (HBRUSH)COLOR_BACKGROUND + 1;
    wc.lpszMenuName  = NULL;
    wc.lpszClassName = APP_NAME;
    wc.hIconSm = (HICON)LoadImage(NULL,
                           MAKEINTRESOURCE(IDI_APPLICATION),
                           IMAGE_ICON,
                           0,
                           0,
                           LR_DEFAULTSIZE | LR_SHARED);
    if(!RegisterClassEx(&wc)) return -1;

    // クライアント領域のサイズから、ウィンドウそのものの幅と高さを計算
    UINT windowWidth  = WIDTH  + GetSystemMetrics(SM_CXFIXEDFRAME) * 2;
    UINT windowHeight = HEIGHT + GetSystemMetrics(SM_CYFIXEDFRAME) * 2
                               + GetSystemMetrics(SM_CYCAPTION);
    // ウィンドウ生成
    HWND hWnd = CreateWindow(APP_NAME,
                    TEXT("はじめてのDelaunay分割"),
                    WS_OVERLAPPEDWINDOW^WS_THICKFRAME^WS_MAXIMIZEBOX,
                    CW_USEDEFAULT,
                    CW_USEDEFAULT,
                    windowWidth,
                    windowHeight,
                    NULL,
                    NULL,
                    hInstance,
                    NULL);
    if(!hWnd) return -1;

    // ウィンドウ表示
    ShowWindow(hWnd, nCmdShow);
    UpdateWindow(hWnd);

    // メッセージループ
    MSG msg;
    while(TRUE)
    {
        if(PeekMessage(&msg, NULL, 0U, 0U, PM_REMOVE))
        {
            if(msg.message == WM_QUIT) break;
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }
    CloseHandle(hMutex);

    return (int)msg.wParam;
}

テスト用の GUI プログラムの方は、思いっきり Windows 依存です。

ただ、Delaunay 分割のアルゴリズムは完全に分離してあるので、他の環境で使いたい場合は Delaunay2d.h だけを持っていって、表示部はがんばって自分で作れば OK です。

肝心の使い方は以下の通りです（Main.cpp のコールバック関数内に実際のサンプルがあります）。

頂点リスト（std::list オブジェクト）を宣言します。
先ほどの頂点リストに、適当な頂点オブジェクトを格納します。
三角形リスト（std::list オブジェクト）を宣言します。
Tercel::Delaunay2d::getDelaunayTriangles() 関数を呼びます。

このとき、getDelaunayTriangles() の第一引数に頂点リストの『参照』を、第2引数には三角形リストの『ポインタ』を渡します。

一見すると奇妙な仕様ですが、関数内で中身が変更されるオブジェクトに関してはポインタ渡し、そうでないものは参照渡しで構文を使い分けているためです。

特に関数を呼び出す側（ここでは Main.cpp ）から見ると、参照渡しは値渡しのように見えるため、関数の中でオブジェクトが変更されないような印象になります。

Delaunay 分割された結果は三角形リストに格納されます。三角形（Tercel::Triangle 構造体）は、メンバに頂点（Tercel::Vertex オブジェクト）のポインタを保持します。ここでは、最初に作った頂点リストの要素を参照する形になります。

ですので、頂点リストの要素を削除した場合、三角形の頂点は不正なポインタになってしまう事に注意する必要があります。

Processing : P3D モードのヘンなクセ

2011-11-05T17:55:00.005+09:00

現行バージョンの Processing には、ソフトウェアレンダリングで 3D を描画できる『P3D モード』が備わっています。

これは、3次元的な情報を視覚化する際にたいへん便利なのですが、実はちょっとした癖もあるため、うまく使うにはそれなりのコツが必要となってきます。

というわけで今日は、P3D モードで遭遇しがちな不具合と、その回避策を紹介したいと思います。

【地面のグリッド表示がおかしくなる現象】

以下のように、地面をグリッドで描画したい場合があったとしましょう。

素直にコードを書くとこんな感じになります（展開してご覧ください）。

void setup() {
  size(800, 600, P3D);
}

void draw() {
  // カメラを適当に設定
  camera(0, -100, 500, 0, 0, 0, 0, 1, 0);
  
  background(255);

  /* ====================== */
  /* 地面の描画（てきとう） */
  /* ====================== */
  stroke(100);
  noFill();  // ←塗りつぶさない
  
  final int step = 20;  
  for(int i = -width; i < width; i += step) {
    beginShape(QUAD_STRIP);
    for(int j = -width; j <= width; j += step) {
      vertex(i, 0, j);
      vertex(i + step, 0, j);
    }
    endShape();
  }
}

上記のプログラムは、いっけん正しく動きそうですが、いざ実行してみるとこんなふうにぶっ飛んだ結果になってしまいます。

なぜこんな事になったのでしょうか。

この現象は、描画すべき線分の端点がカメラ視点の裏側に回り込んでしまった場合に発生します。原因は、Processing 側が 3 次元座標を射影変換する際に、いくつかの例外処理を省略してしまっているためと考えられます。

これを回避するための比較的簡単な方法は、noFill() の代わりに fill(0x00FFFFFF) を使う事です。下記のコードならば、意図した結果が得られるはずです。

void setup() {
  size(800, 600, P3D);
}

void draw() {
  // カメラを適当に設定
  camera(0, -100, 500, 0, 0, 0, 0, 1, 0);
  
  background(255);

  /* ====================== */
  /* 地面の描画（てきとう） */
  /* ====================== */
  stroke(100);
  fill(0x00FFFFFF);  // ←変更
  
  final int step = 20;  
  for(int i = -width; i < width; i += step) {
    beginShape(QUAD_STRIP);
    for(int j = -width; j <= width; j += step) {
      vertex(i, 0, j);
      vertex(i + step, 0, j);
    }
    endShape();
  }
}

基本的に、『塗りつぶし無しの 3D 図形』と『塗りつぶしつきの 3D 図形』では異なるレンダリングパイプラインを通ります。ちなみに、陰影付けやオクルージョンの判定などを必要とする事から、一般的には後者の方が高度な処理となります。

上記のプログラムでは、fill(0x00FFFFFF) によって塗りつぶし色を「完全に透明な白」に指定する事によって、塗りつぶし無しと等価の見た目で、なおかつ問題となっていた端点の座標飛びに対処しています。

【半透明テクスチャのレンダリング順序】

Processing では、透過 png 画像をテクスチャとして使用する事ができます。

これを使えば、以下のように半透明の物体を表現する事ができそうですね。

半透明のテクスチャを貼った四角形を並べるコードは以下のように書けますが、これを実行すると何かがおかしくなります。

void setup() {
  size(800, 600, P3D);
  noStroke();
}

void draw() {
  background(0);
  camera(500, -500, 700, 0, 0, 0, 0, 1, 0);
  
  final int size  = 200;   // 板のサイズ

  final int depth = 1000;  // 板を並べる際の奥行き量
  final int step  =  100;  // 板と板との間隔
  
  colorMode(HSB, depth / step, 100, 100);
  for(int z = depth / 2; z >= -depth / 2; z -= step) {
    
    int colour = color((z + depth/2) / step, 100, 100);  // 色を計算
    PImage tex = createTranslucentTex(colour);           // 指定した色でテクスチャ作成

    // 板の描画
    beginShape(QUADS);
    texture(tex);
    vertex(-size, -size, z, 0,         0);
    vertex(-size,  size, z, 0,         tex.height);
    vertex( size,  size, z, tex.width, tex.height);
    vertex( size, -size, z, tex.width, 0);
    endShape(QUADS);
  }
  colorMode(RGB, 0xFF);

}

// 指定した色で半透明のテクスチャを作成する
// 引数： 色
PImage createTranslucentTex(int colour) {
  PImage tex = createImage(128, 128, ARGB);
  
  for(int y = 0; y < tex.height; y++) {
    int alphaValue = 0xFF * y / tex.height;  // アルファ値
    for(int x = 0; x < tex.width; x++) {
      int index = y * tex.width + x;
      tex.pixels[index] = alphaValue << 24 | 0x00FFFFFF & colour;
    }
  }
  return tex;
}

何がおかしいか分かりますでしょうか。

奥の方の板が、なぜか手前の板を覆い隠しているように見えます（ちなみに、テクスチャが不透明であれば、このような問題は発生しません）。

Processing では、通常『Z バッファリング』と呼ばれる方式を用いてレンダリング処理を省力化しています。それはよいのですが、この方式には、『半透明の物体が必ずしも正しく描画されない』という大きな問題点があります。

これを解決するための手段として、『画家のアルゴリズム（Z ソート法）』というものがあります。アルゴリズムとはいえさして大仰なものではなく、カメラ奥の物体から手前の物体にかけて重ね塗りをするように描画していくという非常に直感的な仕組みです。

というわけで、描画方式をデフォルトのZ バッファリングから画家のアルゴリズムに変更するには、hint() メソッドに ENABLE_DEPTH_SORT を与えればよいそうです。

画家のアルゴリズムを有効にしたコードを以下の通りです。

void setup() {
  size(800, 600, P3D);
  noStroke();
  hint(ENABLE_DEPTH_SORT); // ←追加
}

void draw() {
  background(0);
  camera(500, -500, 700, 0, 0, 0, 0, 1, 0);
  lights();
  final int size  = 200;   // 板のサイズ

  final int depth = 1000;  // 板を並べる際の奥行き量
  final int step  =  100;  // 板と板との間隔

  colorMode(HSB, depth / step, 100, 100);
  for(int z = depth / 2; z >= - depth / 2; z -= step) {
    
    int colour = color((z + depth/2) / step, 100, 100);  // 色を計算
    PImage tex = createTranslucentTex(colour);           // 指定した色でテクスチャ作成

    // 板の描画
    beginShape(QUADS);
    texture(tex);
    vertex(-size, -size, z, 0,         0);
    vertex(-size,  size, z, 0,         tex.height);
    vertex( size,  size, z, tex.width, tex.height);
    vertex( size, -size, z, tex.width, 0);
    endShape(QUADS);
  }
  colorMode(RGB, 0xFF);
}

// 指定した色で半透明のテクスチャを作成する
// 引数： 色
PImage createTranslucentTex(int colour) {
  PImage tex = createImage(128, 128, ARGB);
  
  for(int y = 0; y < tex.height; y++) {
    int alphaValue = 0xFF * y / tex.height;  // アルファ値
    for(int x = 0; x < tex.width; x++) {
      int index = y * tex.width + x;
      tex.pixels[index] = alphaValue << 24 | 0x00FFFFFF & colour;
    }
  }
  return tex;
}

これで、恐らく期待した結果が得られるようになります（『恐らく』と書いたのは、リファレンスに “the algorithm is not yet perfect.” とあるためです）。

【あとがき・おまけ】

実は、オーロラにも半透明のポリゴンが大量に用いられています。ですが、こちらは画家のアルゴリズムを用いていないため、カメラ視点の位置によってはおかしな事になります。

『画家のアルゴリズム』の存在自体は知っていましたが、それを有効にするための方法が分からなかったので、なるべく違和感を払拭するため、ポリゴンの描画順を少しいじったり、視点を移動できる範囲を強く制限したりと涙ぐましい工夫が凝らされています。

というわけで、あまり大した事ではありませんが、知っていると少し幸せになれるかも知れない Tips でした。めでたし、めでたし。

Tercel::Diary

ProcessingでXファイルを解析（スタティックメッシュ篇その2）

ProcessingでXファイルを解析（スタティックメッシュ篇その1）

おもちゃであそぶ

Processingでスキンメッシュアニメーション

Processingで魔法使いの組分け帽子を作る

ProcessingでAI

ProcessingとArduinoをくっつける

ProcessingとTwitterをくっつける

Processingできらきら3次元パーティクル

【予告】 Processingで講演します（？）

Arduinoはじめます

C++と仲直り計画： Cの構造体をスマートポインタで扱う （同時上映： ラムダ式もほんのちょっとかじる）

Processing.jsで作った「ひこうき」を高速化してみた

Processingでてきとう3次元バンプマッピング

三次元・直線・最小自乗法（vs. ぷで）

素晴らしき日々へ

相対参照を適当に考慮したHTMLファイルの移動

冬休みのてきとうアプリケーションづくり

明けましておめでとうございます

Javaでできるだけ安全にデータを直列化するためのてきとう備忘録

近況報告

これなら受かる！ 『ネスペ22 （β版）』

ネットワークスペシャリスト試験に合格しました

Processingでかんたん3Dアニメーション

C++ で 3 次元 Delaunay 分割

続： C++ で Delaunay 分割

C++ で Delaunay 分割（ただし2次元）

Processing : P3D モードのヘンなクセ

C++と仲直り計画： Cの構造体をスマートポインタで扱う（同時上映：ラムダ式もほんのちょっとかじる）

これなら受かる！『ネスペ22 （β版）』