JP7224675B2 - 構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、幾何学的特徴量算出プログラム - Google Patents

Description

本発明は、構造物の形状に関する幾何学的特徴量を算出する構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、特徴量算出プログラムに関するものである。

一般的な計算機援用工学（ＣＡＥ）では、構造物における形状の寸法は、その境界からの距離を計測することにより求められている。

具体的には、特許文献１に示すように、形状の画像データから、画像処理を用いて形状の境界（輪郭）を抽出し、当該境界からの最短距離を形状の寸法としている。

しかしながら、従来の画像処理では、境界を抽出するために近似曲線等の算出が必要であったり、画像データがＳＴＬデータの場合のＳＴＬ要素間の距離の算出が必要であったり、その他、幾何学的情報を得るためには種々の演算処理を行う必要がある。

特開２０１６－１６６８７７号公報

そこで本発明は、上記問題点を解決するためになされたものであり、画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することをその主たる課題としたものである。

すなわち、本発明に係る構造物の幾何学的特徴量算出装置は、構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付ける画像データ受付部と、線形偏微分方程式を用いた仮想的な物理モデルにより、前記画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部と、前記関連値データから前記構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部とを備えることを特徴とする。

また、本発明に係る構造物の幾何学的特徴量算出プログラムは、構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付ける形状データ受付部と、線形偏微分方程式を用いた仮想的な物理モデルにより、前記画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部と、前記関連値データから前記構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部と、としての機能をコンピュータに備えさせることを特徴とする。

このような本発明によれば、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから構造物の幾何学的特徴量に関連する値を示す関連値データを算出し、当該関連値データから構造物の幾何学的特徴量を示す特徴量データを算出することができる。これにより、画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することができる。

関連値データの具体例としては、前記線形偏微分方程式の解を示すデータであることが考えられる。

また、幾何学的特徴量としては、前記構造物の厚さ（幅）、骨格、法線ベクトル又は曲率であることが考えられる。

本発明を用いたシステムとしては、人の身体部位をスキャンして得られる医療用画像データから前記身体部位の幾何学的特徴量を算出する医療用画像分析システム、ハードウエア製品の構造を示す画像データから前記ハードウエア製品の構造の幾何学的特徴量を算出するリバースエンジニアリング支援システム、又は、構造最適化手法により演算装置が演算した最適構造の画像データから前記最適構造の幾何学的特徴量を算出する構造最適化システムなどが考えられる。

構造最適化システムの具体的な実施の態様としては、前記演算装置は、前記幾何学的特徴量算出装置が算出した前記構造物の幾何学的特徴量が所定の幾何学的制約条件を満たすように、トポロジー最適化を行うものであることが考えられる。

積層造形法により造形することのできる構造物を算出するためには、前記演算装置は、前記幾何学的特徴量算出装置が算出した前記構造物の幾何学的特徴量が積層造形法における幾何学的制約条件を満たすように、トポロジー最適化を行うものであることが望ましい。

積層造形法に特有の幾何学的制約条件としては、前記演算装置は、前記積層造形法の造形方向と造形対象の形状の接ベクトルとのなす角を積層傾斜角θ_ａとしたときに、当該積層傾斜角θ_ａを所定角度以下とする幾何学的制約条件を用いるものであることが望ましい。

このように構成した本発明によれば、画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することができる。

本実施形態の構造物の幾何学的特徴量算出装置の機器構成図である。 同実施形態の構造物の幾何学的特徴量算出装置の機能構成図である。 １次元モデルを示す図である。 １次元モデルにおける解析結果を示す図である。 ２次元モデル（リング形状）を示す図である。 ２次元モデル（リング形状）における偏微分方程式の解を示す図である。 ２次元モデル（リング形状）における厚さの出力値を示す図である。 ２次元モデル（リング形状）におけるスケルトンの出力値を示す図である。 ２次元モデル（リング形状）における形状法線ベクトルの出力値を示す図である。 （Ａ）２次元モデル（交差形状）を示す図、（Ｂ）厚さの出力値を示す図、（Ｃ）スケルトンの出力値を示す図である。 （Ａ）２次元モデル（種々の形状）を示す図、（Ｂ）厚さの出力値を示す図、（Ｃ）形状法線ベクトルの出力値を示す図、（Ｄ）スケルトンの出力値を示す図である。 （Ａ）３次元モデルを示す図、（Ｂ）厚さの出力値を示す図である。 （ａ）オーバーハング形状及び積層傾斜角θ_ａを示す図、及び（ｂ）オーバーハング形状における幾何学的特異点を示す図である。 幾何学的特異点の評価方法を説明するための図である。 数値解析に用いた（ａ）モデルを示す図、及び（ｂ）力学的条件を示す図である。 オーバーハング制約を考慮しなかった場合における正則化係数τの各設定値の最適構造を示す図である。 オーバーハング制約を考慮した場合における制御パラメータβの各設定値の最適構造を示す図である。 オーバーハング制約を考慮した場合における各正則化係数τの各設定値の最適構造を示す図である。 オーバーハング制約を考慮した場合における最小オーバーハング角θの各設定値の最適構造を示す図である。 熱の数値解析におけるモデルを示す図である。 （ａ）幾何学的制約を考慮しない場合の最適構造を示す図、及び（ｂ）幾何学的制約を考慮した場合の最適構造を示す図である。

次に、本発明の一実施形態について図面を参照して説明する。なお、図１は本実施形態の構造物の幾何学的特徴量算出装置１００の機器構成図であり、図２は構造物の幾何学的特徴量算出装置１００の機能構成図である。

＜装置構成＞
本実施形態に係る構造物の幾何学的特徴量算出装置１００は、構造物における形状の空間分布の情報のみを用いて、線形偏微分方程式の解から形状の特徴量を解析するものであり、図１に示すように、ＣＰＵ１０１に加えて揮発メモリやＨＤＤ等の記憶装置１０２を備え、さらにマウスやキーボード等の入力手段１０３、演算結果を出力するためのディスプレイやプリンタ等からなる出力手段１０５を接続するための入出力インターフェイス１０４等を有した汎用乃至は専用のコンピュータである。

そして、所定のプログラムを前記記憶装置１０２にインストールし、そのプログラムに基づいてＣＰＵ１０１や周辺機器を協働させることにより、この構造物の特徴量算出装置１００は、図２の機能構成図に示すように、画像データ受付部１、関連値データ算出部２、特徴量データ算出部３、特徴量データ出力部４等としての機能を発揮する。

以下、各部１～４について説明する。

画像データ受付部１は、構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付けるものである。ここで、形状の空間部分を示す画像データとは、空間上における形状のあるなしの分布を示すデータである。また、２次元又は三次元の画像データには、形状のピクセルデータ、ボクセルデータ、例えばＳＴＬ形式のポリゴンデータ等が含まれる。なお、画像データは、例えばユーザが入力手段１０３を用いることにより入力される。

関連値データ算出部２は、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから構造物の幾何学的特徴量に関連する値を示す関連値データを算出するものである。

線形偏微分方程式を用いた仮想的な物理モデルは、以下の式により示される。なお、参照領域Ωは、物体領域Ω_Ｒと空洞領域Ω_Ｒ＼Ωからなる。

ここで、χ_ｉ（∈Ｈ^１（Ω_Ｒ））は、ｉ番目の仮想的な状態変数（偏微分方程式の解）である。ｅ_ｉは、ユークリッド空間Ｒ^ｄにおける標準基底である。ａ^～（チルダ付きのａ）（＞０）は、無次元の拡散係数である。αは、吸収係数である。

また、特性関数１_Ω（∈Ｌ^∞（Ω_Ｒ））は、以下の式により示される。

拡散係数ａ^～として、ａ^～＝ａｈ_０ ^２を満たすように設定している。なお、ｈ_０（＞０）は、代表長さである。また、吸収係数αは、４／ａとしている。なお、吸収係数αは、空洞領域における状態変数χ_ｉをゼロにするために大きな値に設定することが望ましい。

以上から、関連値データ算出部２で用いる仮想的な物理モデルのパラメータは、ａのみとなる。

特徴量データ算出部３は、関連値データ算出部２により関連値データから形状の幾何学的特徴量を示す特徴量データを算出するものである。具体的には、上記の仮想的な物理モデルにより得られた解（χ_ｉ（ｘ））を用いて、形状の幾何学的特徴量を算出する。

特徴量データ算出部３は、上記の解（χ_ｉ（ｘ））を用いて、次式の厚さ関数を定義する。

また、特徴量データ算出部３は、上記の解（χ_ｉ（ｘ））を用いて、次式の形状法線ベクトル関数（形状の外向き法線ベクトルｎを与える関数）を定義する。なお、特徴量データ算出部３は、得られた法線ベクトルｎから、形状の曲率を算出することもできる。

さらに、特徴量データ算出部３は、上記の解（χ_ｉ（ｘ））を用いて、以下により、スケルトン関数を定義する。

まず、形状の幾何学的特徴を表現する形状特徴テンソルＳ^＊を次式により定義する。

ここで、Ｓ^＊の固有値をλ（ｉ）、固有ベクトルをｘ_ｓとする。このとき、ｉはλ（ｉ）≧λ（ｉ＋１）となるように決定する。Ｓ^＊の正規固有ベクトルを用いて基底変換した状態変数をχ^～（チルダ付きのχ）とする。
そして、特徴量データ算出部３は、これを用いて、次式のスケルトン関数を定義する。

ここで、δはディラックのデルタ関数である。

特徴量データ算出部３により得られた上記の特徴量データは、特徴量データ出力部４によってディスプレイに表示される。

＜１次元モデルの解析解＞
図３に示す１次元の材料分布（１次元モデル）において、幾何学的特徴量の算出の妥当性について検討した。

この１次元モデルでは、ｘ軸上において、ｘ＝ｐからｘ＝ｐ＋ｈまで材料が存在している場合を考える。

この１次元モデルの支配方程式（偏微分方程式）及び境界条件は以下となる。

上記の方程式系の解は、以下の式により示される。

ここで、λは√（α／ａ^～）である。

ここで、厚さ関数ｈ（ｘ）は、以下の式により示される。

定数ｃ_ｉは、上記の境界条件、及び、状態変数χとｘ＝ｐ及びｘ＝ｐ＋ｈにおける法線ベクトルの流束に関する連続条件に基づいて定められる。

そして、物体領域が境界ｘ＝Ｋ及びｘ＝Ｌから十分に遠いとすると、

したがって、以下となる。

したがって、厚さ関数ｈ（ｘ）に偏微分方程式の解を代入することによって、１次元モデルの厚さ（長さ）が求まることが分かる。

図４にケース１、ケース２、ケース３として、ｈ、ｐ、ａをそれぞれ変化させた場合の解析結果を示している。

ケース１は、ｈを１．０、１．５、２．０、２．５とした場合であり、ケース２は、ｐを０．３５、０．４、０．４５、０．５とした場合であり、ケース３は、ａを５、１０、２０、５０とした場合である。なお、各ケースにおけるその他のパラメータは以下の表のとおりである。

＜２次元モデル（リング形状）の数値解析結果＞
図５に示す２次元の材料分布（２次元モデル）において幾何学的特徴量の算出の妥当性について検討した。

この２次元モデルは、２次元平面においてリング形状をなすものであり、グレー部分が物体領域（Ω_１～Ω_５）を示し、白い領域が空洞領域を示している。そして、仮想的な物理モデルにおいて、ｈ_０＝０．１とし、ａ＝０．２とした。また、三角形要素を用いて離散化して有限要素法を用いて解析した。

これにより得られた状態変数χ_１、χ_２を図６に示す。また、この状態変数χ_１、χ_２を厚さ関数ｈ（ｘ）に代入することにより、図７に示す厚さ分布データが生成された。

物体領域Ω_１、Ω_２、Ω_３は、互いに異なる厚さを有するリング形状である。各物体領域は、一定の厚さを有している。物体領域Ω_１の厚さは０．１であり、物体領域Ω_２の厚さは０．１２であり、物体領域Ω_３の厚さは０．１４である。したがって、厚さ関数ｈ（ｘ）は、各物体領域において一定値となるはずである。図７に示すとおり、厚さ関数ｈ（ｘ）は、各物体領域において一定値となっていることが分かる。

次に、物体領域Ω_１、Ω_４に着目する。これらの物体領域は、互いに同じ形状であるが、位置が互いに異なるものである。つまり、物体領域Ω_４は、物体領域Ω_１を平行移動させたものである。

図７から分かるように、物体領域Ω_１とΩ_４との厚さ関数ｈ（ｘ）の値は同一であり、同じ形状の物体領域であれば、平行移動させても同じ値を出力することが分かる。

次に、物体領域Ω_１、Ω_５に着目する。これらの物体領域は、互いに同じ形状であるが、位置及び角度が互いに異なるものである。つまり、物体領域Ω_５は、物体領域Ω_１を回転させて平行移動させたものである。

図７から分かるように、物体領域Ω_１とΩ_５との厚さ関数ｈ（ｘ）の値は同一であり、同じ形状の物体領域であれば、回転させて平行移動させても同じ値を出力することが分かる。

また、幾何学的特徴量の算出において厚さ関数とは別の関数（スケルトン関数）を用いることにより、図８に示すように、形状のスケルトン（骨組）を算出することができる。ここで、スケルトン（骨組）は、形状の中心線である。

さらに、幾何学的特徴量の算出において厚さ関数とは別の関数（形状法線ベクトル関数）を用いることにより、図９に示すように、形状の法線ベクトルを算出することができる。

＜２次元モデル（交差形状）の数値解析結果＞
図１０（Ａ）に示す２次元の材料分布（２次元モデル）において幾何学的特徴量の算出の妥当性について検討した。

この２次元モデルは、２次元平面において交差形状をなすものであり、グレー部分が物体領域を示し白い領域が空洞領域を示している。そして、仮想的な物理モデルにおいて、ｈ_０＝０．１とし、ａ＝０．２とした。また、三角形要素を用いて離散化して有限要素法を用いて解析した。

この２次元モデルにおいて得られた厚さ関数ｈ（ｘ）を図１０（Ｂ）に示す。
また、特徴量算出において厚さ関数とは別の関数（スケルトン関数）を用いることにより、図１０（Ｃ）に示すように、形状のスケルトン（骨組）を算出することができる。

その他、図１１（Ａ）に示す２次元の材料分布（２次元モデル）における幾何学的特徴量（厚さ、形状法線ベクトル、スケルトン）を、図１１（Ｂ）～（Ｄ）に示す。

＜３次元モデル（板形状）の数値解析結果＞
図１２（Ａ）に示す３次元の材料分布（３次元モデル）において幾何学的特徴量の算出の妥当性について検討した。

この３次元モデルは、３次元空間において板形状をなすものであり、グレー部分が物体領域を示し、白い領域が空洞領域を示している。そして、仮想的な物理モデルにおいて、ｈ_０＝０．１とし、ａ＝０．２とした。また、テトラ要素を用いて離散化して有限要素法を用いて解析した。

この物理モデルにおいて得られた厚さ関数ｈ（ｘ）を図１２（Ｂ）に示しており、本実施形態により３次元の材料分布においても特徴量（厚さ）を算出できることが分かる。

＜２．本実施形態の効果＞
このように構成した幾何学的特徴量算出装置１００によれば、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから構造物の幾何学的特徴量に関連する値を示す関連値データを算出し、当該関連値データから構造物の幾何学的特徴量を示す特徴量データを算出することができる。これにより、画像データにおける構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することができる。

＜３．本発明の変形実施形態＞
なお、本発明は前記実施形態に限られるものではない。

例えば、前記実施形態の幾何学的特徴量算出装置１００を用いて、医療用画像分析システムを構成しても良い。この場合、画像データ受付部１は、人の身体部位をスキャンして得られる医療用画像データ（例えばＣＴスキャンデータ）を受け付ける。そして、特徴量データ算出部３は、身体部位の幾何学的特徴量を示す特徴量データを算出する。

また、前記実施形態の幾何学的特徴量算出装置１００を用いて、リバースエンジニアリング支援システムを構成しても良い。この場合、画像データ受付部１は、ハードウエア製品の構造を示す画像データを受け付ける。そして、特徴量データ算出部３は、ハードウエア製品の構造の幾何学的特徴量を示す特徴量データを算出する。

さらに、前記実施形態の幾何学的特徴量算出装置１００を用いて、構造最適化システムを構成してもよい。この場合、画像データ受付部１は、構造最適化手法により演算装置が演算した構造物の画像データを受け付ける。ここで、構造最適化手法には、寸法最適化、形状最適化又はトポロジー最適化が考えられる。そして、特徴量データ算出部３は、構造物の幾何学的特徴量を示す特徴量データを算出する。

この構造最適化システムであれば、最適化の過程において形状の寸法を定量的に知ることができるので、構造最適化手法の制約条件として用いることができる。その結果、最大形状寸法、最小形状寸法、又は形状寸法一定等の所定の寸法制約下で構造最適化を行うことができる。構造最適化への適用における最大の特長は、最適化の過程において、必ずしも制約を満たす必要がないので、初期設計解の依存性が低く、高い設計自由度を維持しながら、構造最適化が可能となる。

＜４．構造最適化システム＞
次に、構造最適化システムの一例として、３Ｄプリンタ等に用いられる積層造形法により製造が可能な構造物を創成する構造最適化システムを説明する。

この構造最適化システムは、レベルセット法を用いてトポロジー最適化を行うものである。なお、構造最適化システムは、前記実施形態の幾何学的特徴量算出装置１００が組み込まれたものであり、ＣＰＵに加えて揮発メモリやＨＤＤ等の記憶部を有する演算装置を備え、さらにマウスやキーボード等の入力手段、演算結果を出力するためのディスプレイやプリンタ等からなる出力手段を接続するための入出力インターフェイス等を有した汎用乃至は専用のコンピュータである

＜トポロジー最適化の考え方＞
構造最適化問題とは、ある評価関数Ｊが最小（あるいは最大）となる構造物の形状Ω_１を求める問題である。評価関数は、剛性、揚力、熱伝導特性などの物理的な特性を含む場合が多い。そのため、評価対象とする物理現象を表現する支配方程式を満足することを前提として、最適な形状を求める。支配方程式は、最適化問題における制約条件となり、基本的な構造最適化問題は次のように定式化される。

ここで、ｕは支配方程式の解として得られる状態変数、ｊ（ｕ）は所望の目的関数を表現する関数、Ｕ_ａｄは、許容されるΩ_１がなす集合を表す。評価関数Ｊに含まれる支配方程式の解ｕと形状Ω_１の関係は、陽な形式の関数で表現されない場合が多く、支配方程式の数値解析解として与えられる。次に、構造最適化問題（数 ）に対して、トポロジー最適化の適用について考える。まず、設計の対象とする構造物が配置しうる領域Ω_Ｄ⊂Ｒ^ｎを導入する。ただし、ｎは対象とする空間次元を表す。ここで、領域Ω_Ｄは最適設計解の探索過程において変化しないため、固定設計領域と呼ばれる。固定設計領域は、構造により占められた領域Ω_１（以下、構造領域）と構造に占められていない領域Ω_０（以下、空洞領域）により構成されており、それらを識別する特性関数χを次式により定義する。

ただし、構造領域Ω_１と空洞領域Ω_０の境界は構造領域に含めるものとする。このとき、特性関数χを用いて、トポロジー最適化問題は次のように定式化される。

ただし，Ｆ_ａｄは許容される形状の特性関数のなす集合を表す。トポロジー最適化問題では、設計変数が特性関数χであるから、最適構造設計問題が材料の最適分布問題として表現される特徴を持つ。その結果、設計変更の過程において、構造領域Ω_１のトポロジーが変更される可能性を持つ。

＜トポロジー最適化問題の緩和と正則化＞
トポロジー最適化問題は、最適解が存在しないため、不良設定問題である(Allaire, 2002)。従って、何らかの手段を用いて良設定問題に置き換える手続を必要とする。設計空間の緩和に基づく方法として、均質化設計法(Bendsoe and Kikuchi, 1988)が代表的である。また、広く普及している密度法(Bendsoe and Sigmund, 1999)は、均質化設計法の考え方に基づいて簡便な問題に置き換えた方法である。これらの方法では、構造領域と空洞領域の中間状態を許容するため、連続な密度関数を用いて最適構造が表現される。他方、正則化を用いた方法としては、レベルセット法に基づく方法(山田他,2009; Yamada et al., 2010)がある。この方法では，レベルセット関数と呼ばれるスカラー関数の等位面により構造領域の境界面を表現し、レベルセット関数の変更により構造領域の形状変更を表現する方法である。また、レベルセット関数に滑らかさを適切に確保することで、トポロジー最適化問題の正則化を行う。ここでは、外形形状の幾何学的制約について扱うため、レベルセット法に基づく方法を用いる。

＜レベルセット法に基づく方法＞
レベルセット法に基づく方法では、次式に示すスカラー関数φ（ｘ）を導入し、形状表現を行う。

レベルセット関数を用いて特性関数χを再定義すれば、次式となる。

レベルセット関数により定義された特性関数χを用いてトポロジー最適化問題を定式化する。

次に、レベルセット関数の更新方法について述べる。基本的な考え方は。最適化の過程におけるレベルセット関数φ（ｘ）の変更を、仮想的な時間ｔにおける時間変化とすることである。レベルセット関数φ（ｘ）は、値域制約（－１≦φ（ｘ）≦１）を満たしながら、次式に示す反応拡散方程式に基づいて更新される。

ただし、Ｋ∈Ｒ_＋は比例定数、Ｃ∈Ｒ_＋は正規化定数、τ∈Ｒ_＋は正則化定数、Ｊ’はトポロジー導関数(Amstutz et al.,2006)である。反応拡散方程式における拡散項は、レベルセット関数φ（ｘ）の滑らかさを保証するだけではなく、拡散項の相対的な寄与度を与える正則化係数τの調節により、得られる最適構造の幾何学的複雑さの定性的な制御ができる特徴も持つ。さらには、正則化係数の寄与を方向毎に異なる大きさに設定すれば，幾何学的複雑さに指向性を持たせることが可能となる。その結果、フライス加工等に要求される等断面形状制約を容易に考慮できる(Yamada et al., 2010)。また、比例定数Ｋと正規化定数Ｃは、レベルセット関数の値域制約が１程度であるから、同程度の大きさに設定するパラメータである。

＜幾何学的特徴量に対する制約を考慮したトポロジー最適化＞
＜幾何学的特徴量に対する偏微分方程式系＞
幾何学的特徴量を抽出する対象形状Ω_１を包含する領域Ωにおいて、次に示す偏微分方程式系を考える。

ここで，ｓ_ｉ∈Ｈ^１（Ω）はｘ_ｉ方向に対する状態変数、ｅ_ｉは標準基底、ａ∈Ｒ_＋は拡散係数である。なお，拡散係数ａは相対的な影響が小さくなるように、無次元化拡散係数ａ^～：＝ａ／Ｌ^２を十分小さな値ａ^～≪１に仮定する。ただし、Ｌは代表長さである。ここで、上記の偏微分方程式の力学的意味を議論するために、状態変数ｓ_ｉから構成されるベクトル場ｓ＝［ｓ_１ｓ_２・・・ｓ_ｎ］^Ｔ及びその試験関数∀ｓ^～∈Ｈ_０ ^１（Ω）^ｎを導入すると、次式の等価な弱形式が得られる。

ベクトル場ｓの支配方程式は、縦弾性係数がａ、ポアソン比が０である等方性線形弾性体で占められた領域Ωに対して、周囲∂Ωを完全拘束、領域Ω_０に大きさが１の等方性ばねを配置、対象形状の境界∂Ω_１から外向きに大きさ１の界面力を作用させた場合の弾性方程式である。従って、ベクトル場ｓは、このような力学的条件における変位場と等価である。直感的には、ポアソン比が０であるから、ベクトル場ｓの向きは。境界∂Ω_１において、その法線ベクトルの向きと一致すると推察される。従って、対象形状Ω_１の法線ベクトルを近似的に表現する関数ｆｎ（｛ｓ_ｉ｝_{１≦ｉ≪ｎ}；ε_ｓｎ）は次式で与えられる。

ただし，ε_ｓｎ∈Ｒ_＋は、ｓ＝０の場合においてもベクトル関数ｆｎが定義できるように導入された微小な正の値を持つパラメータである。なお、ｓ＝０となるのは、境界面∂Ω_１から十分に遠い場所であるから、微小パラメータε_ｓｎの設定値がベクトル関数ｆｎの意味に実質的に影響を与えないことを注記する。形状の境界面∂Ω_１における厳密な法線ベクトルは、正則な境界面上のみで定義されるが、ベクトル関数ｆｎは、接ベクトルの勾配が不連続となる幾何学的特異点を含めて定義できる。従って、ベクトル関数ｆｎは、正則な境界面近傍においては、法線ベクトルに漸近する関数であり、本来の定義より広い意味で法線ベクトル定義し、領域Ω全域で値を持つ。証明には、偏微分方程式の解の振る舞いに関する厳密な議論が必要であり、今後の検討課題ではあるが、数値解析に基づく検証の結果、高い精度で近似されることが確認されている(Yamada, 2019b)。また、偏微分方程式における拡散項の特徴から、半径がＯ（√ａ）の範囲で状態変数ｓ_ｉが平均化された振る舞いとなる。つまり、ベクトル関数ｆｎは、∂Ω_１の近傍において、半径Ｏ（√ａ）の大きさで平均化された法線ベクトルであると言える。

次に、状態変数ｓ_ｉの分布に関する基本的な特徴を議論するために、２次元円形形状に対する偏微分方程式の解析解を示す。円の半径をＲとし、円の中心を原点とするとき、解析解は、極座標表示（ｒ，θ）を用いて、次式となる。

ただし、Ｋ_ｍ（ｘ）は第二種変形ベッセル関数である。上式（数２３）に示すように、形状Ω_１の内部においては、原点を通る線形関数、外部では急激かつ単調に減小しながら零に収束する関数である。したがって、ベクトル場ｓの大きさの絶対値は。境界∂Ω_１に近づくにつれて単調に増加する。また、ベクトル場ｓの向きは、境界∂Ω_１において、境界の法線方向と厳密に一致している。領域Ω_１の内部と外部においても、最も近い境界面の法線方向と一致している。

＜オーバーハング制約に対する幾何学的制約条件＞
次に、積層造形法におけるオーバーハング形状に関する幾何学的制約条件を明確化する。最初に、図１３（ａ）に示すように、造形方向と造形対象の形状（以下、造形形状）の接ベクトルとのなす角を積層傾斜角θａと定義する。積層傾斜角θａが一定値を超えると、造形中の自重の作用により、造形物の表面精度が著しく悪化する問題がある。例えば、図１３（ａ）の場合、破線で示した造形形状は、積層造形法による製造が困難であると言える。このように、下部に十分な支持がないために製造が困難な部分形状をオーバーハング形状、その部分形状の水平面からの傾斜角θをオーバーハング角と呼ぶ。この問題を回避するために、オーバーハング形状の下部に、サポート材と呼ばれる仮の造形物を配置する対処法が用いられる。しかしながら、サポート材の配置に基づく対処法は、製造工程やコストの増加を招くだけではなく、サポート材の除去が可能な単純な造形形状に限定される。

そこで、サポート材なしに積層造形を可能とする造形形状の幾何学的条件を明確化する。先行研究の多くは、形状の法線ベクトルに基づいて幾何学的制約条件を定式化している。法線ベクトルに基づく考え方を図１３（ｂ）に示す造形形状に適用すれば、破線で示した境界が制約違反と判断される。しかしながら、点Ｐ_ｓ ^１＋と点Ｐ_ｓ ^２＋及びその上部は、明らかに積層造形が不可能な形状である。

本実施形態では、造形が成立しない領域の評価及び平均曲率に基づいた幾何学的特異点の評価に基づく方法を提案する。具体的には、図１４（ａ）に示す円形形状において，積層方向をｘ_２軸方向と設定した場合、境界ｌ_ｃを評価するのではなく、境界ｌ_ｃが円弧となる扇形の領域は、制約違反により成立する領域（以下、制約違反領域）として定義する。この領域は、図１４（ｂ）に示す左側への積層造形ができない領域（以下、負側違反領域）と、図１４（ｃ）に示す右側への積層造形ができない領域（以下、正側違反領域）のいずれにも該当する領域である。なお、３次元の場合は、法線ベクトルを水平面へ投影した場合の向き及び鉛直方向が基底ベクトルとなる２次元平面を考え、同様の議論を考えれば良い。

次に、図１３（ｂ）の点Ｐ_ｓ ^１＋と点Ｐ_ｓ ^２＋のような幾何学的特異点の抽出法について考える。ここで、点Ｐ_ｓ ^１＋と点Ｐ_ｓ ^２＋のような凸型形状の幾何学的特異点では、平均曲率は正の無限大であり，点Ｐ_ｓ ^１－と点Ｐ_ｓ ^２－のような凹型形状の幾何学的特異点では平均曲率は負の無限大であることに注目する。平均曲率は法線ベクトルの発散により与えられるから、幾何学的特異点の近傍領域は、ベクトル関数ｆｎの発散の絶対値が大きい領域であり、その正負により凸型形状もしくは凹型形状であるかを識別できる。なお、前述のように、ベクトル関数ｆｎはＯ（√ａ）程度の範囲で平均化されるため、幾何学的特異点近傍においても、ベクトル関数ｆｎの発散は有界な値を持つことを注記する。

以上の議論から、積層造形におけるオーバーハング形状に対する幾何学的制約条件は、制約違反領域及び制約違反領域近傍の凸型形状の幾何学的特異点が空集合となるように定式化すれば良い。

＜幾何学的制約条件の定式化＞
幾何学的特徴量に対する偏微分方程式系及び上記の議論に基づき、２次元問題における幾何学的制約条件を定式化する。単純化のために、積層方向がｘ２の正の方向となるように座標系を選ぶこととする。このとき、積層限界方向を与える単位ベクトルξ_１及びξ_２は、幾何学的関係から次のように与えられる。

負側違反領域は、ベクトルξ_１と拡張した法線ベクトルｆ_ｎの内積が正となる対象形状Ω_１の内部であるから、ベクトル場ｓは対象形状Ω_１近傍のみで非零であることに注目すれば、領域Ω全体における負側違反領域近傍は、ベクトルξ_１とベクトル場ｓの内積が正の領域として与えられる。同様に、正側違反領域は、ベクトルξ_２とベクトル場ｓの内積が正となる領域として与えられる。従って、制約違反領域に基づいた幾何学的制約条件は次式により与えられる。

ただし、Ｒはランプ関数である。なお、ランプ関数の代わりにヘビサイド関数を用いた定式化も可能であるが、ヘビサイド関数の導関数はデルタ関数であるから、最適化問題の多峰性を強めてしまう問題を持つことを注記する。さらに、ベクトル場ｓはＯ（１／√ａ）の大きさを持つことに注目して適切に正規化を行うと、幾何学的制約条件は次のよう定式化される。

ここで，Ｔ（ｘ；γ）は次式である。

パラメータγを十分に大きく設定した場合、上式の被積分関数値は、制約違反領域の内部において一定（Ｔ＝１）となる。なお、空間分布を正規化するＴは、次に示す四つの条件を満たす必要がある。

次に、幾何学的特異点について考える。前述のように、拡張された法線ベクトルの発散の絶対値が大きい制約違反領域近傍が幾何学的特異点の近傍である。さらに、その符号が正の領域は凸型形状であるから、制約汎関数値Ｇ＾_ａｍの被積分関数を過大に評価すれば良い。同様に符号が負の領域は凹型形状であるから、制約汎関数値Ｇ＾_ａｍの被積分関数を過小評価すれば良い。ただし、負の値をとることを許容してはならない。従って、幾何学特異点を考慮した制約汎関数は次式となる。

ここで、β∈Ｒ_＋は、幾何学的特異点の排除度合いを与えるパラメータである。パラメータβの設定値が大きい場合、制約違反に近い形状、かつ、局所的に曲率が大きい場所を排除する。そのため，許容されるオーバーハング角の制約を違反する部分構造を許容する。従って、幾何学的特異点に対する制御パラメータβは、下向き凸型形状の幾何学的特異点を排除しながら、適切な曲率を維持できる範囲に設定する必要がある。パラメータβの設定値の具体的な影響については、後述の数値解析例で示す。

＜最適化問題の定式化＞
次に、弾性テンソルがＤの等方性線形弾性体で構成される構造領域Ω_１について考える。構造領域Ω_１は、全境界∂Ωの一部の境界Γ_ｕで変位拘束、別の一部の境界Γ_ｔで表面力ｔ∈Ｌ^２（Γ_ｔ）^２が負荷されているとする。このとき、変位場をｕ∈Ｈ_Γｕ ^１（Ω_１）^２とし、平面ひずみを仮定された二次元問題における体積制約及びオーバーハング制約を考慮した平均コンプライアンス最小化問題は次のように定式化される。

ただし、ｈは領域の厚み、Ｇ_ｍａｘ∈Ｒ_＋は許容される体積の上限値、ε（ｕ）：＝１／２×（∇ｕ＋∇ｕ^Ｔ）はひずみであり、Ｈ_Γｕ ^１は次に示す関数空間である。

また，弾性テンソルＤは、縦弾性係数Ｅとポアソン比ｖを用いて次式で与えられる。

ただし、δ_ｉｊはクロネッカーのデルタである。

＜数値実装法＞
＜最適化アルゴリズム＞
構造最適化システムの演算装置は、次に示すアルゴリズムに従って最適化を行う。
（ステップ１）
レベルセット関数の適当な初期値φ_０を設定する。
（ステップ２）
有限要素法を用いて、変位場及びベクトル場ｓを数値解析する。
（ステップ３）
目的汎関数及び制約汎関数を計算し、収束していれば最適化を終了する。
（ステップ４）
有限要素法を用いてベクトル場ｓに対する随伴場を数値解析し、随伴変数法により設計感度（トポロジー導関数）を求める。
（ステップ５）
反応拡散方程式を用いてレベルセット関数φを更新し、（ステップ２）へ戻る。
以上の手順により最適構造の候補を示すレベルセット関数値を得る。

上記（ステップ１）のレベルセット関数は、演算装置のレベルセット関数データ格納部に格納されている。レベルセット関数は、初期構造が設定された設計領域の各部が、構造を形成し、物体により占められた物体領域（物体相）、空洞を形成する空洞領域（空洞相）、又はそれら領域の境界であるかを示すものであり、物体領域を表す値及び空洞領域を表す値の間の所定値が、物体領域及び空洞領域の境界を表す。なお、レベルセット関数データは、例えばユーザが入力手段を用いることにより入力される。

また、設計領域を示す設計領域データは、演算装置の設計領域データ格納部に格納されている。なお、設計領域（当該設計領域を要素分割する構造格子（メッシュ）情報を含む。）を示す設計領域データは、例えばユーザが入力手段を用いることにより入力される。

さらに、前記設計領域（固定設計領域）の境界条件を示す境界条件データは、演算装置の境界条件データ格納部に格納されている。具体的な境界条件としては、例えば設計領域の拘束条件、初期構造に作用する荷重等の外力（表面力）などである。なお、境界条件データは、例えばユーザが入力手段を用いることにより入力される。

その他、最適構造を求めるときに用いる体積制約に必要な制約の条件値及び変位場の解析に必要なヤング率やポアソン比の値等の材料定数である解析データは、演算装置の解析データ格納部に格納されている。なお、解析データは、例えばユーザが入力手段を用いることにより入力される。

そして、（ステップ２）～（ステップ５）は、演算装置のレベルセット関数更新部により実行される。レベルセット関数更新部は、所定の制約条件下において、剛性、固有振動数等の構造物の性能を目標値に近づけるようにレベルセット関数を更新して、物体領域と空洞領域との境界を移動させるとともに、レベルセット関数の更新に伴う物体領域内のトポロジー変化（形態変化）を許容して物体領域内に新たな空洞領域（穴）を形成し、その新たな空洞領域と物体領域との境界を移動させるものである。

＜オイラー座標系に基づく変位場の近似解法＞
構造領域Ω_１はオイラー座標系により表現されており、最適化の過程における反復計算ごとに、有限要素の生成が必要となる。計算コストの観点から、Ersatz material approach（Allaire et al., 2004）による近似解法を適用する。すなわち、空洞領域Ω_０は相対的に小さい弾性テンソルを持つ構造材料と仮定し、その界面近傍では滑らかに分布する材料特性を持つと仮定する。具体的には、有限要素解析において、次に示す拡張した弾性テンソルＤ^～を用いて、変位場に対する支配方程式を固定設計領域Ω_Ｄに対して拡張する。

ただし、ε_ｕ∈Ｒ_＋は、近似ヘビサイド関数Ｈ^～の遷移幅を与える微小パラメータ、ｄ∈Ｒ_＋は、構造領域の縦弾性係数に対する空洞領域の相対値である。

＜ベクトル場ｓに対する特性関数の拡張＞
ベクトル場ｓは固定設計領域Ω_Ｄを包含する領域Ωにおいて解析する必要があるため、領域Ωへ広げた拡張平滑化特性関数Ψ∈Ｈ^１（Ω）^ｎを次式により定義する。

ただし、ε_Ψ∈Ｒ_＋は、平滑化度合いを与える微小パラメータ、ε_ｗ∈Ｒ_＋は、特性関数の遷移幅を与える微小パラメータである。

ランプ関数Ｒ（ｘ）はｘ＝０において強い意味で微分が定義できないため、微小パラメータε_ｓｒ∈Ｒ_＋を導入し、ランプ関数Ｒに漸近する双曲線を近似ランプ関数Ｒ^～として定義する。

設計感度を導出する感度解析の段階において、近似ランプ関数を用いて数値解析を行う。

＜数値解析例＞
数値解析に用いたモデルは、図１５（ａ）に示すように、６×１の矩形領域を固定設計領域とし、この下辺と同一直線上に積層プラットフォームを設定し、積層方向は鉛直上向きとする。このとき、縦弾性係数が１、ポアソン比が０．３の平面ひずみを仮定された等方性線形弾性体で占められた形状Ω_１の体積制約付き剛性最大化問題について考える。

力学的条件は、図１５（ｂ）に示すように、下辺両端の境界Γ_ｕにおいてローラー支持し、上辺中央下向きに大きさｔ＝１の表面力を作用させる、左右対称であるから、中央に対称条件を与えることで、右半分の領域のみを解析する。ベクトル場ｓを定義する領域Ωは、固定設計領域を包含するように、左辺、右辺及び上端側に幅０．２を持つ矩形領域分を広げた領域とした。下辺側に超過した制約違反領域を評価する必要がないため、領域を広げる必要はないことを注記する。また、許容される体積の上限値は固定設計領域の５０％する。固定設計領域は、１５０×５０の構造格子に対して、格子中心周りに対称４分割した３角形要素を用いて有限要素分割し、変位場、幾何学状態場及びレベルセット関数場を１次有限要素により離散化する。反応拡散方程式及び幾何学的特徴量に対する偏微分方程式の無次元化に用いる代表長さＬはＬ＝１とする。また、無次元化拡散係数ａ^～は幾何学的特徴量の評価をする際の平滑化半径の大きさを与えるから、幾何学的複雑さを与える正則化係数τと比較して十分に小さく設定する必要がある。ここでは、ａ^～＝０．０３τと設定する．オーバーハング制約汎関数Ｇ_ａｍの正規化パラメータγはγ＝２．５、ランプ関数の平滑化パラメータε_ｓｒは、適切に滑らかさを保持できるように、ξ・ｓに対してはε_ｓｒ＝３、曲率評価（１＋β×ｄｉｖｆｎ）に対してはε_ｓｒ＝１×１０^－９とする。変位場の近似解法に用いるパラメータは、ｄ＝１×１０^－３、ε_ｕ＝０．５とする．幾何学的特徴量に対する偏微分方程式の近似解法に用いるパラメータは、ε_Ψ＝１×１０^－４Ｌ^２、ε_ｗ＝０．９５とする。

最初に、オーバーハング制約を考慮しなかった場合の最適構造を図１６に示す。図中、黒色で占められた領域が構造領域Ω_１、灰色で絞められた領域が空洞領域Ω_０を表す。図１６に示すように、幾何学的正脚条件を考慮しない場合、積層フラットフォームと平行に近い部分構造を持つ最適構造が得られる。また、正則化係数τは構造の幾何学的複雑さを制御する効果を持つものの、正則化係数τの調整を行っても、積層造形が困難な部分構造を排除できないことを確認できる。

＜幾何学的特異点に対する制御パラメータβの影響＞
幾何学的特異点に対する制御パラメータβが最適構造に与える影響について検討する。許容される最小のオーバーハング角をθ＝４５°に設定し、幾何学特異点に対する制御パラメータβをβ＝３×１０^－２、２×１０^－２、１×１０^－２、８×１０^－３、３×１０^－３、５×１０^－４、１×１０^－４、０に設定した。得られた最適構造を図１７に示す。β＝０の場合、幾何学的特異点に対する考慮がないため、図１７（ｈ）に示すように、凸型の幾何学的特異点を持つ部分構造が得られている。制御パラメータβの設定値が大きい場合、図１７（ａ）に示すように，許容されるオーバーハング角の制約を違反する部分構造の領域が広いことが確認できる。制御パラメータβが小さくなるにつれて、許容される曲率の範囲が広くなるため、凹型の幾何学的特異点に近づいていくことも確認できる。従って、制御パラメータβを適切な範囲に設定すれば、オーバーハング制約を満たさない部分領域が非常に小さく、積層造形が可能な最適構造が得られることがわかった。また、幾何学的特異点を考慮しなければ、所望のオーバーハング制約を満たす最適構造が得られないこともわかった。
なお、以降の数値解析例においては、β＝６×１０^－３に設定する。

＜正則化係数τの影響＞
次に、異なる正則化係数τを設定し、最適化を図った。許容される最小のオーバーハング角はθ＝４５°とする。このとき、得られた最適構造を図１８に示す。図１８に示すように、いずれの構造も、所望のオーバーハング制約を満たす構造が得られていることが確認できる。また、正則化係数τの設定値の大きさにより、最適構造の幾何学的な複雑さが異なることを確認できる。従って、オーバーハング制約を考慮した場合においても、最適構造の幾何学的複雑さを定性的に制御可能であることがわかった。

＜最小許容オーバーハング角の影響＞
最後に、許容される最小オーバーハング角の設定値が最適構造に与える影響について検討する。許容される最小のオーバーハング角として、４つの設定値（θ１５°、３０°、４５°、６０°）について最適化を図った。正則化係数は、τ＝３×１０^－３及びτ＝４×１０^－４に設定し、得られた最適構造を図１９に示す。図１９に示すように、いずれの結果も、オーバーハング制約を考慮した結果が得られ、オーバーハング角の設定値毎に異なる最適構造が得られた。

上記の構造最適化システムは、オーバーハング形状に関する幾何学的制約条件を加味したトポロジー最適化を行うものであったが、幾何学的制約条件としては、オーバーハング形状に限られない。例えば、図２０に示すモデルを用いた体積制約付き剛性最大化問題において、構造物の厚さ（幅）を所定値以下とする幾何学的制約条件を加味したトポロジー最適化を行うようにしてもよい。その結果を図２１に示す。

その他、本発明は前記実施形態に限られず、その趣旨を逸脱しない範囲で種々の変形が可能であるのは言うまでもない。

１００・・・構造物の幾何学的特徴量算出装置
１ ・・・画像データ受付部
２ ・・・関連値データ算出部
３ ・・・特徴量データ算出部
４ ・・・特徴量データ出力部

Claims (10)

1. 構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付ける画像データ受付部と、
   線形偏微分方程式を用いた仮想的な物理モデルにより、前記画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部と、
   前記関連値データから前記構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部とを備える、構造物の幾何学的特徴量算出装置。
2. 前記関連値データは、前記線形偏微分方程式の解を示すデータである、請求項１記載の構造物の幾何学的特徴量算出装置。
3. 前記幾何学的特徴量は、前記構造物の厚さ、骨格、法線ベクトル又は曲率である、請求項１記載の構造物の幾何学的特徴量算出装置。
4. 請求項１記載の構造物の幾何学的特徴量算出装置を用いて構成され、人の身体部位をスキャンして得られる医療用画像データから前記身体部位の幾何学的特徴量を算出する医療用画像分析システム。
5. 請求項１記載の構造物の幾何学的特徴量算出装置を用いて構成され、ハードウエア製品の構造を示す画像データから前記ハードウエア製品の構造の幾何学的特徴量を算出するリバースエンジニアリング支援システム。
6. 請求項１記載の構造物の幾何学的特徴量算出装置を用いて構成され、構造最適化手法により演算装置が演算した構造物の画像データから前記構造物の幾何学的特徴量を算出する構造最適化システム。
7. 前記演算装置は、前記幾何学的特徴量算出装置が算出した前記構造物の幾何学的特徴量が所定の幾何学的制約条件を満たすように、トポロジー最適化を行うものである、請求項６記載の構造最適化システム。
8. 前記演算装置は、前記幾何学的特徴量算出装置が算出した前記構造物の幾何学的特徴量が積層造形法における幾何学的制約条件を満たすように、トポロジー最適化を行うものである、請求項６記載の構造最適化システム。
9. 前記演算装置は、前記積層造形法の造形方向と造形対象の形状の接ベクトルとのなす角を積層傾斜角θ_ａとしたときに、当該積層傾斜角θ_ａを所定角度以下とする幾何学的制約条件を用いるものである、請求項８記載の構造最適化システム。
10. 構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付ける形状データ受付部と、
    線形偏微分方程式を用いた仮想的な物理モデルにより、前記画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部と、
    前記関連値データから前記構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部と、としての機能をコンピュータに備えさせる、構造物の幾何学的特徴量算出プログラム。

Images