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

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

Info

Publication number
WO2019216221A1
WO2019216221A1 PCT/JP2019/017419 JP2019017419W WO2019216221A1 WO 2019216221 A1 WO2019216221 A1 WO 2019216221A1 JP 2019017419 W JP2019017419 W JP 2019017419W WO 2019216221 A1 WO2019216221 A1 WO 2019216221A1
Authority
WO
WIPO (PCT)
Prior art keywords
geometric feature
geometric
feature amount
shape
value
Prior art date
Application number
PCT/JP2019/017419
Other languages
English (en)
French (fr)
Inventor
崇恭 山田
Original Assignee
国立大学法人京都大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立大学法人京都大学 filed Critical 国立大学法人京都大学
Priority to JP2020518248A priority Critical patent/JP7224675B2/ja
Publication of WO2019216221A1 publication Critical patent/WO2019216221A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes

Definitions

  • the present invention relates to a structure geometric feature amount calculation device, a medical image analysis system, a reverse engineering support system, a structure optimization system, and a feature amount calculation program for calculating a geometric feature amount related to the shape of a structure. It is about.
  • CAE computer aided engineering
  • a shape boundary (contour) is extracted from shape image data using image processing, and the shortest distance from the boundary is defined as the shape dimension.
  • the present invention has been made to solve the above problems, and its main problem is to calculate the geometric feature amount of the structure without extracting the boundary of the shape of the structure from the image data. It is what.
  • the geometric feature amount calculation apparatus for a structure uses an image data receiving unit that receives two-dimensional or three-dimensional image data indicating a spatial distribution of a shape in the structure, and a linear partial differential equation.
  • a related value data calculating unit for calculating related value data indicating a value related to a geometric feature of the structure from the image data by a virtual physical model; and a geometric value of the structure from the related value data.
  • a feature amount data calculation unit for calculating feature amount data indicating the target feature amount.
  • the geometric feature amount calculation program for a structure uses a shape data receiving unit that receives two-dimensional or three-dimensional image data indicating a spatial distribution of the shape in the structure, and a linear partial differential equation.
  • a related value data calculating unit for calculating related value data indicating a value related to a geometric feature of the structure from the image data by a virtual physical model; and a geometric value of the structure from the related value data.
  • the computer is provided with a function as a feature amount data calculation unit that calculates feature amount data indicating a characteristic feature amount.
  • the related value data indicating the value related to the geometric feature amount of the structure is calculated from the image data by the virtual physical model using the linear partial differential equation, and the related value Feature quantity data indicating the geometric feature quantity of the structure can be calculated from the data.
  • the geometric feature amount of the structure can be calculated without extracting the boundary of the shape of the structure from the image data.
  • the related value data it may be data indicating a solution of the linear partial differential equation.
  • the geometric feature amount may be the thickness (width), skeleton, normal vector, or curvature of the structure.
  • a medical image analysis system for calculating a geometric feature amount of a body part from medical image data obtained by scanning a human body part, an image showing a structure of a hardware product
  • the reverse engineering support system that calculates the geometric feature value of the hardware product structure from the data, or the geometric feature value of the optimal structure from the image data of the optimal structure calculated by the calculation device using the structure optimization method
  • a structural optimization system to be calculated can be considered.
  • the arithmetic device may be configured such that the geometric feature amount of the structure calculated by the geometric feature amount calculation device satisfies a predetermined geometric constraint condition.
  • topology optimization is performed.
  • the arithmetic device is configured such that the geometric feature amount of the structure calculated by the geometric feature amount calculating device is a geometric feature in the additive manufacturing method. It is desirable to perform topology optimization so as to satisfy the physical constraints.
  • the geometric feature amount of the structure can be calculated without extracting the boundary of the shape of the structure from the image data.
  • FIG. A) The figure which shows a three-dimensional model, (B) The figure which shows the output value of thickness.
  • overhang shape shows the (a) overhang shape and stacked tilt angle theta a. It is a figure for demonstrating the evaluation method of a geometric singularity. It is a figure which shows the (a) model used for numerical analysis, and (b) is a figure which shows mechanical conditions. It is a figure which shows the optimal structure of each setting value of the regularization coefficient (tau) when an overhang restriction is not considered. It is a figure which shows the optimal structure of each setting value of the control parameter (beta) in the case of considering overhang restrictions. It is a figure which shows the optimal structure of each setting value of each regularization coefficient (tau) in the case of considering overhang restrictions.
  • FIG. 1 is a device configuration diagram of a structural geometric feature amount calculating apparatus 100 according to the present embodiment
  • FIG. 2 is a functional configuration diagram of the structural geometric feature amount calculating apparatus 100.
  • the geometric feature value calculation device 100 for a structure analyzes a feature value of a shape from a solution of a linear partial differential equation using only information on the spatial distribution of the shape in the structure.
  • a storage device 102 such as a volatile memory or an HDD is provided.
  • an input unit 103 such as a mouse or a keyboard, and an output unit 105 including a display or a printer for outputting calculation results are provided.
  • It is a general purpose or dedicated computer having an input / output interface 104 or the like for connection.
  • the structure feature amount calculating apparatus 100 is configured as shown in the functional configuration diagram of FIG. In addition, it functions as the image data receiving unit 1, the related value data calculating unit 2, the feature amount data calculating unit 3, the feature amount data output unit 4, and the like.
  • the image data receiving unit 1 receives two-dimensional or three-dimensional image data indicating the spatial distribution of the shape in the structure.
  • the image data indicating the space portion of the shape is data indicating a distribution with or without a shape in the space.
  • the two-dimensional or three-dimensional image data includes shape pixel data, voxel data, such as polygon data in the STL format. Note that the image data is input by the user using the input unit 103, for example.
  • the related value data calculation unit 2 calculates related value data indicating a value related to the geometric feature amount of the structure from the image data by a virtual physical model using a linear partial differential equation.
  • the reference region ⁇ includes an object region ⁇ R and a cavity region ⁇ R ⁇ ⁇ .
  • ⁇ i ( ⁇ H 1 ( ⁇ R )) is the i-th virtual state variable (solution of the partial differential equation).
  • e i is a standard basis in the Euclidean space R d .
  • a ⁇ (a with tilde) (> 0) is the diffusion coefficient of dimensionless.
  • ⁇ is an absorption coefficient.
  • the absorption coefficient ⁇ is 4 / a.
  • the absorption coefficient ⁇ is preferably set to a large value in order to make the state variable ⁇ i in the cavity region zero.
  • the parameter of the virtual physical model used in the related value data calculation unit 2 is only a.
  • the feature value data calculation unit 3 is used by the related value data calculation unit 2 to calculate feature value data indicating the geometric feature value of the shape from the related value data. Specifically, the geometric feature quantity of the shape is calculated using the solution ( ⁇ i (x)) obtained by the virtual physical model.
  • the feature amount data calculation unit 3 defines a thickness function of the following equation using the above solution ( ⁇ i (x)).
  • the feature amount data calculation unit 3 defines a shape normal vector function (a function that gives an outward normal vector n of the shape) using the above solution ( ⁇ i (x)).
  • the feature amount data calculation unit 3 can also calculate the curvature of the shape from the obtained normal vector n.
  • the feature amount data calculation unit 3 uses the above solution ( ⁇ i (x)) to define a skeleton function as follows.
  • the shape feature tensor S * expressing the geometric feature of the shape is defined by the following equation.
  • i is determined so that ⁇ (i) ⁇ ⁇ (i + 1).
  • a state variable subjected to basis transformation using the normal eigenvector of S * is defined as ⁇ ⁇ ( ⁇ with tilde).
  • the feature-value data calculation part 3 defines the skeleton function of following Formula using this.
  • ⁇ is a Dirac delta function
  • the feature amount data obtained by the feature amount data calculation unit 3 is displayed on the display by the feature amount data output unit 4.
  • ⁇ is the ⁇ ( ⁇ / a ⁇ ).
  • the thickness function h (x) is expressed by the following equation.
  • the thickness (length) of the one-dimensional model can be obtained by substituting the solution of the partial differential equation into the thickness function h (x).
  • FIG. 4 shows the analysis results when h, p, and a are changed as Case 1, Case 2, and Case 3, respectively.
  • Case 1 is a case where h is 1.0, 1.5, 2.0, and 2.5
  • Case 2 is a case where p is 0.35, 0.4, 0.45, and 0.5
  • Case 3 is a case where a is set to 5, 10, 20, 50. The other parameters in each case are as shown in the table below.
  • This two-dimensional model has a ring shape on a two-dimensional plane, and a gray portion indicates an object region ( ⁇ 1 to ⁇ 5 ), and a white region indicates a hollow region.
  • FIG. 6 shows the state variables ⁇ 1 and ⁇ 2 obtained thereby. Further, the thickness distribution data shown in FIG. 7 is generated by substituting the state variables ⁇ 1 and ⁇ 2 into the thickness function h (x).
  • the object regions ⁇ 1 , ⁇ 2 , and ⁇ 3 are ring shapes having different thicknesses. Each object region has a certain thickness.
  • the thickness of the object region ⁇ 1 is 0.1
  • the thickness of the object region ⁇ 2 is 0.12
  • the thickness of the object region ⁇ 3 is 0.14. Therefore, the thickness function h (x) should be a constant value in each object region. As shown in FIG. 7, it can be seen that the thickness function h (x) has a constant value in each object region.
  • the object regions ⁇ 1 and ⁇ 4 These object areas have the same shape but different positions. That is, the object region ⁇ 4 is obtained by translating the object region ⁇ 1 .
  • the values of the thickness functions h (x) of the object regions ⁇ 1 and ⁇ 4 are the same, and if the object regions have the same shape, the same value is output even if they are translated. I understand.
  • the object regions ⁇ 1 and ⁇ 5 These object regions have the same shape but different positions and angles. That is, the object region ⁇ 5 is obtained by rotating and moving the object region ⁇ 1 .
  • the values of the thickness function h (x) of the object regions ⁇ 1 and ⁇ 5 are the same, and if the object region has the same shape, the same value can be obtained even if it is rotated and translated. You can see that it outputs.
  • a shape skeleton (frame) can be calculated.
  • the skeleton is the center line of the shape.
  • shape normal vector function shape normal vector function
  • This two-dimensional model has an intersecting shape in a two-dimensional plane, and a gray portion indicates an object region and a white region indicates a hollow region.
  • the thickness function h (x) obtained in this two-dimensional model is shown in FIG. Further, by using a function (skeleton function) different from the thickness function in the feature amount calculation, as shown in FIG. 10C, a shape skeleton (frame) can be calculated.
  • 11A to 11D show geometric features (thickness, shape normal vector, skeleton) in the two-dimensional material distribution (two-dimensional model) shown in FIG. 11A.
  • the three-dimensional model has a plate shape in a three-dimensional space, and a gray portion indicates an object region and a white region indicates a hollow region.
  • the thickness function h (x) obtained in this physical model is shown in FIG. 12B, and it can be seen that the feature quantity (thickness) can be calculated even in the three-dimensional material distribution according to this embodiment.
  • the related value indicating the value related to the geometric feature amount of the structure from the image data by the virtual physical model using the linear partial differential equation.
  • Data can be calculated, and feature amount data indicating the geometric feature amount of the structure can be calculated from the related value data.
  • the geometric feature amount of the structure can be calculated without extracting the boundary of the shape of the structure in the image data.
  • a medical image analysis system may be configured using the geometric feature value calculation apparatus 100 of the above-described embodiment.
  • the image data receiving unit 1 receives medical image data (for example, CT scan data) obtained by scanning a human body part.
  • the feature-value data calculation part 3 calculates the feature-value data which show the geometric feature-value of a body part.
  • a reverse engineering support system may be configured using the geometric feature value calculation apparatus 100 of the above embodiment.
  • the image data receiving unit 1 receives image data indicating the structure of the hardware product.
  • the feature amount data calculation unit 3 calculates feature amount data indicating the geometric feature amount of the structure of the hardware product.
  • a structure optimization system may be configured using the geometric feature value calculation apparatus 100 of the above embodiment.
  • the image data receiving unit 1 receives the image data of the structure calculated by the calculation device using the structure optimization method.
  • dimensional optimization, shape optimization, or topology optimization can be considered as the structure optimization method.
  • the feature amount data calculation unit 3 calculates feature amount data indicating the geometric feature amount of the structure.
  • This structure optimization system can quantitatively know the size of the shape in the optimization process, and can be used as a constraint condition for the structure optimization method.
  • structural optimization can be performed under predetermined dimensional constraints such as maximum shape size, minimum shape size, or constant shape size.
  • the biggest feature of application to structural optimization is that it is not always necessary to satisfy the constraints in the optimization process, so the initial design solution is less dependent and structural optimization is possible while maintaining high design flexibility. It becomes.
  • the structure optimization system performs topology optimization using the level set method.
  • the structure optimization system incorporates the geometric feature quantity calculation device 100 of the above-described embodiment, and includes an arithmetic device having a storage unit such as a volatile memory and an HDD in addition to a CPU, and a mouse, A general purpose or dedicated computer having an input means such as a keyboard, an input / output interface for connecting output means such as a display or a printer for outputting calculation results, etc.
  • the structure optimization problem is a problem for obtaining the shape ⁇ 1 of the structure in which a certain evaluation function J is minimum (or maximum).
  • the evaluation function often includes physical characteristics such as rigidity, lift, and heat conduction characteristics. Therefore, an optimum shape is obtained on the assumption that the governing equation expressing the physical phenomenon to be evaluated is satisfied.
  • the governing equation becomes a constraint condition in the optimization problem, and the basic structural optimization problem is formulated as follows.
  • u is a state variable obtained as a solution of the governing equation
  • j (u) is a function representing a desired objective function
  • U ad is a set formed by allowed ⁇ 1 .
  • the relationship between the solution u of the governing equation included in the evaluation function J and the shape ⁇ 1 is often not expressed by an explicit function, and is given as a numerical analysis solution of the governing equation.
  • the fixed design region is composed of a region ⁇ 1 occupied by the structure (hereinafter referred to as a structural region) and a region ⁇ 0 not occupied by the structure (hereinafter referred to as a hollow region). It is defined by an expression.
  • the boundary between the structural region ⁇ 1 and the cavity region ⁇ 0 is included in the structural region.
  • the topology optimization problem is formulated as follows.
  • Fad represents a set of characteristic functions of allowable shapes.
  • the design variable is the characteristic function ⁇
  • the optimum structure design problem has a feature expressed as an optimum material distribution problem.
  • the topology of the structural region ⁇ 1 is changed in the process of the design change.
  • the topology optimization problem is a defect setting problem because there is no optimal solution (Allaire, 2002). Therefore, a procedure for replacing the good setting problem with some means is required.
  • the homogenized design method (Bendsoe and Kikuchi, 1988) is a typical method based on relaxation of the design space.
  • the widely used density method (Bendsoe and Sigmund, 1999) is a method that is replaced with a simple problem based on the concept of the homogenization design method. In these methods, an optimum structure is expressed using a continuous density function in order to allow an intermediate state between the structure region and the cavity region.
  • a method using regularization there is a method based on the level set method (Yamada et al., 2009; Yamada et al., 2010).
  • a boundary surface of a structural region is expressed by a coordinate surface of a scalar function called a level set function
  • a shape change of the structural region is expressed by changing the level set function.
  • the topology optimization problem is regularized by appropriately ensuring smoothness in the level set function.
  • a method based on the level set method is used in order to deal with geometric constraints on the outer shape.
  • the change of the level set function ⁇ (x) in the optimization process is a time change at the virtual time t.
  • the level set function ⁇ (x) is updated based on the reaction diffusion equation shown below, while satisfying the range constraint ( ⁇ 1 ⁇ ⁇ (x) ⁇ 1).
  • K ⁇ R + is a proportionality constant
  • C ⁇ R + is a normalization constant
  • ⁇ R + is a regularization constant
  • J ′ is a topological derivative
  • the proportionality constant K and the normalization constant C are parameters that are set to the same magnitude because the range restriction of the level set function is about 1.
  • s i ⁇ H 1 ( ⁇ ) is a state variable in the x i direction
  • e i is a standard basis
  • a ⁇ R + is a diffusion coefficient.
  • the diffusion coefficient a as relative impact is reduced, the non-dimensional diffusion coefficients a ⁇ : assume a / L 2 to a sufficiently small value a ⁇ ⁇ 1.
  • L is a representative length.
  • the vector field s [s 1 s 2 ... S n ] T composed of the state variables s i and its test functions ⁇ s ⁇ ⁇ If H 0 1 ( ⁇ ) n is introduced, an equivalent weak form of the following equation is obtained.
  • the governing equation of the vector field s, longitudinal elastic modulus a, with respect to Poisson's ratio 0 a is isotropic linear elastic body in occupied area Omega, fully constrain the periphery ⁇ Omega, the size in the region Omega 0
  • This is an elastic equation when an isotropic spring of 1 is arranged and an interface force of a magnitude 1 is applied outward from the boundary ⁇ 1 of the target shape. Therefore, the vector field s is equivalent to the displacement field in such a mechanical condition.
  • the direction of the vector field s is.
  • a function fn ( ⁇ s i ⁇ 1 ⁇ i ⁇ n ; ⁇ sn ) that approximately represents the normal vector of the target shape ⁇ 1 is given by the following equation.
  • the s 0 since it is sufficiently far away from the boundary surface ⁇ Omega 1, noted that the setting value of the minute parameter epsilon sn does not substantially affect the meaning of the vector function fn.
  • the exact normal vector at the shape boundary ⁇ 1 is defined only on the regular boundary surface, but the vector function fn is defined including the geometric singularity where the gradient of the tangent vector becomes discontinuous. it can.
  • the vector function fn is a function that asymptotically approaches the normal vector in the vicinity of the regular boundary surface, and the normal vector is defined in a broader sense than the original definition and has a value in the entire region ⁇ .
  • Proof requires a strict discussion on the behavior of the solution of the partial differential equation, and although it is a subject for future study, as a result of verification based on numerical analysis, it has been confirmed that it can be approximated with high accuracy (Yamada , 2019b).
  • the behavior is such that the state variable s i is averaged within a radius of O ( ⁇ a). That is, the vector function fn can be said to be a normal vector averaged with the size of the radius O ( ⁇ a) in the vicinity of ⁇ 1 .
  • K m (x) is a second type modified Bessel function.
  • the shape ⁇ 1 has a linear function passing through the origin and a function that converges to zero while decreasing rapidly and monotonously outside. Therefore, the absolute value of the magnitude of the vector field s. As it approaches the boundary ⁇ 1 , it increases monotonously.
  • the direction of the vector field s is the boundary ⁇ Omega 1, match exactly with the normal direction of the boundary. Also in internal and external areas Omega 1, coincides with the normal direction of the nearest boundary surface.
  • ⁇ Geometric constraints for overhang constraints an angle formed by a modeling direction and a tangent vector of a shape to be modeled (hereinafter, modeling shape) is defined as a stacking inclination angle ⁇ a. If the stacking inclination angle ⁇ a exceeds a certain value, there is a problem that the surface accuracy of the modeled object is remarkably deteriorated due to the action of its own weight during modeling. For example, in the case of FIG. 13A, it can be said that the modeling shape indicated by the broken line is difficult to manufacture by the layered modeling method.
  • an overhang shape a partial shape that is difficult to manufacture due to lack of sufficient support at the lower portion
  • an inclination angle ⁇ from the horizontal plane of the partial shape is called an overhang angle.
  • a countermeasure is used in which a temporary model called a support material is arranged at the lower part of the overhang shape.
  • the handling method based on the arrangement of the support material is not limited to a manufacturing process and cost, but is limited to a simple modeling shape that can remove the support material.
  • constraint violation area an area established by constraint violation (hereinafter referred to as constraint violation area).
  • This area includes an area where additive manufacturing to the left side shown in FIG. 14B (hereinafter referred to as a negative violation area) and an area where additive manufacturing to the right side shown in FIG. 14C cannot be performed (hereinafter referred to as positive violation).
  • This is a region corresponding to any of (regions).
  • a similar discussion may be considered by considering a two-dimensional plane in which the normal vector is projected onto a horizontal plane and the direction and the vertical direction are base vectors.
  • the average curvature is positive infinity, such as the point P s 1 ⁇ and the point P s 2 ⁇ .
  • the average curvature is negative infinity at the geometric singularity of the concave shape. Since the average curvature is given by the divergence of the normal vector, the region near the geometric singularity is a region where the absolute value of the divergence of the vector function fn is large, and whether the convex shape or the concave shape is positive or negative. Can be identified. Note that, as described above, since the vector function fn is averaged in the range of about O ( ⁇ a), the divergence of the vector function fn has a bounded value even near the geometric singularity. To do.
  • the geometric constraint condition for the overhang shape in additive manufacturing may be formulated so that the geometric singularities of the convex shape in the vicinity of the constraint violation region and the constraint violation region become an empty set.
  • the vector field s is non-zero only in the vicinity of the target shape ⁇ 1. If attention is paid, the vicinity of the negative violation region in the entire region ⁇ is given as a region where the inner product of the vector ⁇ 1 and the vector field s is positive. Similarly, the positive violation region is given as a region where the inner product of the vector ⁇ 2 and the vector field s is positive. Accordingly, the geometric constraint condition based on the constraint violation region is given by the following equation.
  • R is a ramp function. It should be noted that although a formulation using a snake side function can be used instead of a ramp function, since the derivative of the snake side function is a delta function, it has a problem of increasing the multimodality of the optimization problem. . Further, if the vector field s is appropriately normalized by paying attention to the magnitude of O (1 / ⁇ a), the geometric constraint condition is formulated as follows.
  • T (x; ⁇ ) is the following equation.
  • the vicinity of the constraint violation area where the absolute value of the divergence of the extended normal vector is large is the vicinity of the geometric singularity. Furthermore, since the area
  • ⁇ R + is a parameter that gives the degree of exclusion of geometric singularities.
  • ⁇ R + is a parameter that gives the degree of exclusion of geometric singularities.
  • h is the thickness of the region
  • G max ⁇ R + is the upper limit of the allowable volume
  • ⁇ (u): 1 ⁇ 2 ⁇ ( ⁇ u + ⁇ u T ) is the strain
  • H ⁇ u 1 is This is the function space shown.
  • the elastic tensor D is given by the following equation using the longitudinal elastic modulus E and Poisson's ratio v.
  • Step 1 Set an appropriate initial value ⁇ 0 of the level set function.
  • Step 2 Using the finite element method, the displacement field and the vector field s are numerically analyzed.
  • Step 3 The objective functional and the constraint functional are calculated, and if they have converged, the optimization ends.
  • Step 4 The adjoint field for the vector field s is numerically analyzed using the finite element method, and the design sensitivity (topological derivative) is obtained by the adjoint variable method.
  • Step 5 Update the level set function ⁇ using the reaction diffusion equation, and return to (Step 2).
  • the level set function value indicating the optimum structure candidate is obtained by the above procedure.
  • the level set function of (Step 1) above is stored in the level set function data storage unit of the arithmetic unit.
  • the level set function is a structure in which each part of the design area in which the initial structure is set forms a structure, the object area occupied by the object (object phase), the cavity area forming the cavity (cavity phase), or the boundary between these areas
  • the predetermined value between the value representing the object region and the value representing the cavity region represents the boundary between the object region and the cavity region.
  • the level set function data is input by the user using an input unit, for example.
  • design area data indicating the design area is stored in the design area data storage unit of the arithmetic device.
  • design area data indicating a design area including structural grid (mesh) information that divides the design area into elements
  • an input unit for example.
  • boundary condition data indicating the boundary condition of the design area is stored in the boundary condition data storage unit of the arithmetic unit.
  • Specific boundary conditions include, for example, constraint conditions in the design region, external force (surface force) such as a load acting on the initial structure, and the like.
  • the boundary condition data is input by the user using an input unit, for example.
  • the analytical data that is the material constants such as the value of the constraint necessary for volume constraint used when obtaining the optimal structure and the value of Young's modulus and Poisson's ratio necessary for the analysis of the displacement field are stored in the analytical data storage unit of the arithmetic unit Stored in The analysis data is input by the user using an input unit, for example.
  • Step 2 to (Step 5) are executed by the level set function update unit of the arithmetic unit.
  • the level set function update unit updates the level set function so that the performance of the structure such as rigidity and natural frequency approaches the target value under a predetermined constraint condition, and moves the boundary between the object region and the cavity region
  • the topology change (morphological change) in the object region with the update of the level set function is allowed, and a new cavity region (hole) is formed in the object region, and the boundary between the new cavity region and the object region Is to move.
  • the structural region ⁇ 1 is expressed by the Euler coordinate system, and it is necessary to generate a finite element for each iterative calculation in the optimization process. From the viewpoint of calculation cost, the approximate solution method by Arbitr material approach (Allaire et al., 2004) is applied. That is, it is assumed that the cavity region ⁇ 0 is a structural material having a relatively small elastic tensor, and has material characteristics that are smoothly distributed in the vicinity of the interface. Specifically, in the finite element analysis using the elastic tensor D ⁇ an extension shown below, to extend the governing equation for the displacement field with respect to the fixed design domain Omega D.
  • epsilon u ⁇ R + is small parameter that gives an approximate Heaviside function H transition width of ⁇
  • d ⁇ R + is the relative value of the cavity region relative to the longitudinal elastic modulus of the structural region.
  • ⁇ ⁇ ⁇ R + is a minute parameter that gives the degree of smoothing
  • ⁇ w ⁇ R + is a minute parameter that gives the transition width of the characteristic function.
  • ⁇ Example of numerical analysis> In the model used for the numerical analysis, as shown in FIG. 15A, a 6 ⁇ 1 rectangular area is set as a fixed design area, a stacking platform is set on the same straight line as the lower side, and the stacking direction is vertically upward. .
  • a volume-constrained stiffness maximization problem of shape ⁇ 1 occupied by an isotropic linear elastic body assuming a plane strain with a longitudinal elastic modulus of 1 and a Poisson's ratio of 0.3.
  • the region ⁇ defining the vector field s is a region obtained by expanding a rectangular region having a width of 0.2 on the left side, the right side, and the upper end side so as to include the fixed design region. Note that there is no need to widen the area because there is no need to evaluate the constraint violation area that exceeded the lower side. Further, the upper limit value of the allowable volume is 50% of the fixed design area.
  • the fixed design area is divided into finite elements using a triangular element divided into four symmetrically around the lattice center of a 150 ⁇ 50 structured lattice, and the displacement field, geometric state field, and level set function field are first-order finite. Discretize by elements.
  • ⁇ sr 1 ⁇ 10 ⁇ 9 .
  • FIG. 16 shows an optimum structure when the overhang constraint is not taken into consideration.
  • the region occupied by black represents the structural region ⁇ 1
  • the region confined in gray represents the cavity region ⁇ 0 .
  • the regularization coefficient ⁇ has an effect of controlling the geometric complexity of the structure, it can be confirmed that even if the regularization coefficient ⁇ is adjusted, a partial structure that is difficult to be layered cannot be excluded.
  • control parameter ⁇ for the geometric singularity The effect of the control parameter ⁇ for the geometric singularity on the optimal structure is discussed.
  • ⁇ Influence of minimum allowable overhang angle> Finally, the influence of the set value of the minimum allowable overhang angle on the optimum structure is examined.
  • the minimum allowable overhang angle was optimized for four set values ( ⁇ 15 °, 30 °, 45 °, 60 °).
  • topology optimization performs topology optimization in consideration of the geometric constraint condition regarding the overhang shape, but the geometric constraint condition is not limited to the overhang shape.
  • topology optimization may be performed in consideration of a geometric constraint condition in which the thickness (width) of the structure is a predetermined value or less. Good. The result is shown in FIG.

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Image Analysis (AREA)

Abstract

画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出する。構造物の２次元又は３次元の画像データを受け付ける画像データ受付部１と、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部２と、関連値データから構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部３とを備える。

Description

構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、幾何学的特徴量算出プログラム
　本発明は、構造物の形状に関する幾何学的特徴量を算出する構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、特徴量算出プログラムに関するものである。
　一般的な計算機援用工学（ＣＡＥ）では、構造物における形状の寸法は、その境界からの距離を計測することにより求められている。
　具体的には、特許文献１に示すように、形状の画像データから、画像処理を用いて形状の境界（輪郭）を抽出し、当該境界からの最短距離を形状の寸法としている。
　しかしながら、従来の画像処理では、境界を抽出するために近似曲線等の算出が必要であったり、画像データがＳＴＬデータの場合のＳＴＬ要素間の距離の算出が必要であったり、その他、幾何学的情報を得るためには種々の演算処理を行う必要がある。
特開２０１６－１６６８７７号公報
　そこで本発明は、上記問題点を解決するためになされたものであり、画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することをその主たる課題としたものである。
　すなわち、本発明に係る構造物の幾何学的特徴量算出装置は、構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付ける画像データ受付部と、線形偏微分方程式を用いた仮想的な物理モデルにより、前記画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部と、前記関連値データから前記構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部とを備えることを特徴とする。
　また、本発明に係る構造物の幾何学的特徴量算出プログラムは、構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付ける形状データ受付部と、線形偏微分方程式を用いた仮想的な物理モデルにより、前記画像データから前記構造物の幾何学的特徴量に関連する値を示す関連値データを算出する関連値データ算出部と、前記関連値データから前記構造物の幾何学的特徴量を示す特徴量データを算出する特徴量データ算出部と、としての機能をコンピュータに備えさせることを特徴とする。
　このような本発明によれば、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから構造物の幾何学的特徴量に関連する値を示す関連値データを算出し、当該関連値データから構造物の幾何学的特徴量を示す特徴量データを算出することができる。これにより、画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することができる。
　関連値データの具体例としては、前記線形偏微分方程式の解を示すデータであることが考えられる。
　また、幾何学的特徴量としては、前記構造物の厚さ（幅）、骨格、法線ベクトル又は曲率であることが考えられる。
　本発明を用いたシステムとしては、人の身体部位をスキャンして得られる医療用画像データから前記身体部位の幾何学的特徴量を算出する医療用画像分析システム、ハードウエア製品の構造を示す画像データから前記ハードウエア製品の構造の幾何学的特徴量を算出するリバースエンジニアリング支援システム、又は、構造最適化手法により演算装置が演算した最適構造の画像データから前記最適構造の幾何学的特徴量を算出する構造最適化システムなどが考えられる。
　構造最適化システムの具体的な実施の態様としては、前記演算装置は、前記幾何学的特徴量算出装置が算出した前記構造物の幾何学的特徴量が所定の幾何学的制約条件を満たすように、トポロジー最適化を行うものであることが考えられる。
　積層造形法により造形することのできる構造物を算出するためには、前記演算装置は、前記幾何学的特徴量算出装置が算出した前記構造物の幾何学的特徴量が積層造形法における幾何学的制約条件を満たすように、トポロジー最適化を行うものであることが望ましい。
　積層造形法に特有の幾何学的制約条件としては、前記演算装置は、前記積層造形法の造形方向と造形対象の形状の接ベクトルとのなす角を積層傾斜角θａとしたときに、当該積層傾斜角θａを所定角度以下とする幾何学的制約条件を用いるものであることが望ましい。
　このように構成した本発明によれば、画像データにおいて構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することができる。
本実施形態の構造物の幾何学的特徴量算出装置の機器構成図である。 同実施形態の構造物の幾何学的特徴量算出装置の機能構成図である。 １次元モデルを示す図である。 １次元モデルにおける解析結果を示す図である。 ２次元モデル（リング形状）を示す図である。 ２次元モデル（リング形状）における偏微分方程式の解を示す図である。 ２次元モデル（リング形状）における厚さの出力値を示す図である。 ２次元モデル（リング形状）におけるスケルトンの出力値を示す図である。 ２次元モデル（リング形状）における形状法線ベクトルの出力値を示す図である。 （Ａ）２次元モデル（交差形状）を示す図、（Ｂ）厚さの出力値を示す図、（Ｃ）スケルトンの出力値を示す図である。 （Ａ）２次元モデル（種々の形状）を示す図、（Ｂ）厚さの出力値を示す図、（Ｃ）形状法線ベクトルの出力値を示す図、（Ｄ）スケルトンの出力値を示す図である。 （Ａ）３次元モデルを示す図、（Ｂ）厚さの出力値を示す図である。 （ａ）オーバーハング形状及び積層傾斜角θａを示す図、及び（ｂ）オーバーハング形状における幾何学的特異点を示す図である。 幾何学的特異点の評価方法を説明するための図である。 数値解析に用いた（ａ）モデルを示す図、及び（ｂ）力学的条件を示す図である。 オーバーハング制約を考慮しなかった場合における正則化係数τの各設定値の最適構造を示す図である。 オーバーハング制約を考慮した場合における制御パラメータβの各設定値の最適構造を示す図である。 オーバーハング制約を考慮した場合における各正則化係数τの各設定値の最適構造を示す図である。 オーバーハング制約を考慮した場合における最小オーバーハング角θの各設定値の最適構造を示す図である。 熱の数値解析におけるモデルを示す図である。 （ａ）幾何学的制約を考慮しない場合の最適構造を示す図、及び（ｂ）幾何学的制約を考慮した場合の最適構造を示す図である。
　次に、本発明の一実施形態について図面を参照して説明する。なお、図１は本実施形態の構造物の幾何学的特徴量算出装置１００の機器構成図であり、図２は構造物の幾何学的特徴量算出装置１００の機能構成図である。
　＜装置構成＞
　本実施形態に係る構造物の幾何学的特徴量算出装置１００は、構造物における形状の空間分布の情報のみを用いて、線形偏微分方程式の解から形状の特徴量を解析するものであり、図１に示すように、ＣＰＵ１０１に加えて揮発メモリやＨＤＤ等の記憶装置１０２を備え、さらにマウスやキーボード等の入力手段１０３、演算結果を出力するためのディスプレイやプリンタ等からなる出力手段１０５を接続するための入出力インターフェイス１０４等を有した汎用乃至は専用のコンピュータである。
　そして、所定のプログラムを前記記憶装置１０２にインストールし、そのプログラムに基づいてＣＰＵ１０１や周辺機器を協働させることにより、この構造物の特徴量算出装置１００は、図２の機能構成図に示すように、画像データ受付部１、関連値データ算出部２、特徴量データ算出部３、特徴量データ出力部４等としての機能を発揮する。
　以下、各部１～４について説明する。
　画像データ受付部１は、構造物における形状の空間分布を示す２次元又は３次元の画像データを受け付けるものである。ここで、形状の空間部分を示す画像データとは、空間上における形状のあるなしの分布を示すデータである。また、２次元又は三次元の画像データには、形状のピクセルデータ、ボクセルデータ、例えばＳＴＬ形式のポリゴンデータ等が含まれる。なお、画像データは、例えばユーザが入力手段１０３を用いることにより入力される。
　関連値データ算出部２は、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから構造物の幾何学的特徴量に関連する値を示す関連値データを算出するものである。
　線形偏微分方程式を用いた仮想的な物理モデルは、以下の式により示される。なお、参照領域Ωは、物体領域ΩＲと空洞領域ΩＲ＼Ωからなる。
Figure JPOXMLDOC01-appb-M000001
　ここで、χｉ（∈Ｈ１（ΩＲ））は、ｉ番目の仮想的な状態変数（偏微分方程式の解）である。ｅｉは、ユークリッド空間Ｒｄにおける標準基底である。ａ～（チルダ付きのａ）（＞０）は、無次元の拡散係数である。αは、吸収係数である。
　また、特性関数１Ω（∈Ｌ∞（ΩＲ））は、以下の式により示される。
Figure JPOXMLDOC01-appb-M000002
　拡散係数ａ～として、ａ～＝ａｈ０ ２を満たすように設定している。なお、ｈ０（＞０）は、代表長さである。また、吸収係数αは、４／ａとしている。なお、吸収係数αは、空洞領域における状態変数χｉをゼロにするために大きな値に設定することが望ましい。
　以上から、関連値データ算出部２で用いる仮想的な物理モデルのパラメータは、ａのみとなる。
　特徴量データ算出部３は、関連値データ算出部２により関連値データから形状の幾何学的特徴量を示す特徴量データを算出するものである。具体的には、上記の仮想的な物理モデルにより得られた解（χｉ（ｘ））を用いて、形状の幾何学的特徴量を算出する。
　特徴量データ算出部３は、上記の解（χｉ（ｘ））を用いて、次式の厚さ関数を定義する。
Figure JPOXMLDOC01-appb-M000003
　また、特徴量データ算出部３は、上記の解（χｉ（ｘ））を用いて、次式の形状法線ベクトル関数（形状の外向き法線ベクトルｎを与える関数）を定義する。なお、特徴量データ算出部３は、得られた法線ベクトルｎから、形状の曲率を算出することもできる。
Figure JPOXMLDOC01-appb-M000004
　さらに、特徴量データ算出部３は、上記の解（χｉ（ｘ））を用いて、以下により、スケルトン関数を定義する。
　まず、形状の幾何学的特徴を表現する形状特徴テンソルＳ＊を次式により定義する。
Figure JPOXMLDOC01-appb-M000005
 　ここで、Ｓ＊の固有値をλ（ｉ）、固有ベクトルをｘｓとする。このとき、ｉはλ（ｉ）≧λ（ｉ＋１）となるように決定する。Ｓ＊の正規固有ベクトルを用いて基底変換した状態変数をχ～（チルダ付きのχ）とする。
　そして、特徴量データ算出部３は、これを用いて、次式のスケルトン関数を定義する。
Figure JPOXMLDOC01-appb-M000006
 　ここで、δはディラックのデルタ関数である。
　特徴量データ算出部３により得られた上記の特徴量データは、特徴量データ出力部４によってディスプレイに表示される。
　＜１次元モデルの解析解＞
　図３に示す１次元の材料分布（１次元モデル）において、幾何学的特徴量の算出の妥当性について検討した。
　この１次元モデルでは、ｘ軸上において、ｘ＝ｐからｘ＝ｐ＋ｈまで材料が存在している場合を考える。
　この１次元モデルの支配方程式（偏微分方程式）及び境界条件は以下となる。
Figure JPOXMLDOC01-appb-M000007
　上記の方程式系の解は、以下の式により示される。
Figure JPOXMLDOC01-appb-M000008
　ここで、λは√（α／ａ～）である。
　ここで、厚さ関数ｈ（ｘ）は、以下の式により示される。
Figure JPOXMLDOC01-appb-M000009
　定数ｃｉは、上記の境界条件、及び、状態変数χとｘ＝ｐ及びｘ＝ｐ＋ｈにおける法線ベクトルの流束に関する連続条件に基づいて定められる。
Figure JPOXMLDOC01-appb-M000010
　そして、物体領域が境界ｘ＝Ｋ及びｘ＝Ｌから十分に遠いとすると、
Figure JPOXMLDOC01-appb-M000011
　したがって、以下となる。
Figure JPOXMLDOC01-appb-M000012
Figure JPOXMLDOC01-appb-M000013
　したがって、厚さ関数ｈ（ｘ）に偏微分方程式の解を代入することによって、１次元モデルの厚さ（長さ）が求まることが分かる。
　図４にケース１、ケース２、ケース３として、ｈ、ｐ、ａをそれぞれ変化させた場合の解析結果を示している。
　ケース１は、ｈを１．０、１．５、２．０、２．５とした場合であり、ケース２は、ｐを０．３５、０．４、０．４５、０．５とした場合であり、ケース３は、ａを５、１０、２０、５０とした場合である。なお、各ケースにおけるその他のパラメータは以下の表のとおりである。
Figure JPOXMLDOC01-appb-T000014
　＜２次元モデル（リング形状）の数値解析結果＞
　図５に示す２次元の材料分布（２次元モデル）において幾何学的特徴量の算出の妥当性について検討した。
　この２次元モデルは、２次元平面においてリング形状をなすものであり、グレー部分が物体領域（Ω１～Ω５）を示し、白い領域が空洞領域を示している。そして、仮想的な物理モデルにおいて、ｈ０＝０．１とし、ａ＝０．２とした。また、三角形要素を用いて離散化して有限要素法を用いて解析した。
　これにより得られた状態変数χ１、χ２を図６に示す。また、この状態変数χ１、χ２を厚さ関数ｈ（ｘ）に代入することにより、図７に示す厚さ分布データが生成された。
　物体領域Ω１、Ω２、Ω３は、互いに異なる厚さを有するリング形状である。各物体領域は、一定の厚さを有している。物体領域Ω１の厚さは０．１であり、物体領域Ω２の厚さは０．１２であり、物体領域Ω３の厚さは０．１４である。したがって、厚さ関数ｈ（ｘ）は、各物体領域において一定値となるはずである。図７に示すとおり、厚さ関数ｈ（ｘ）は、各物体領域において一定値となっていることが分かる。
　次に、物体領域Ω１、Ω４に着目する。これらの物体領域は、互いに同じ形状であるが、位置が互いに異なるものである。つまり、物体領域Ω４は、物体領域Ω１を平行移動させたものである。
　図７から分かるように、物体領域Ω１とΩ４との厚さ関数ｈ（ｘ）の値は同一であり、同じ形状の物体領域であれば、平行移動させても同じ値を出力することが分かる。
　次に、物体領域Ω１、Ω５に着目する。これらの物体領域は、互いに同じ形状であるが、位置及び角度が互いに異なるものである。つまり、物体領域Ω５は、物体領域Ω１を回転させて平行移動させたものである。
　図７から分かるように、物体領域Ω１とΩ５との厚さ関数ｈ（ｘ）の値は同一であり、同じ形状の物体領域であれば、回転させて平行移動させても同じ値を出力することが分かる。
　また、幾何学的特徴量の算出において厚さ関数とは別の関数（スケルトン関数）を用いることにより、図８に示すように、形状のスケルトン（骨組）を算出することができる。ここで、スケルトン（骨組）は、形状の中心線である。
　さらに、幾何学的特徴量の算出において厚さ関数とは別の関数（形状法線ベクトル関数）を用いることにより、図９に示すように、形状の法線ベクトルを算出することができる。
　＜２次元モデル（交差形状）の数値解析結果＞
　図１０（Ａ）に示す２次元の材料分布（２次元モデル）において幾何学的特徴量の算出の妥当性について検討した。
　この２次元モデルは、２次元平面において交差形状をなすものであり、グレー部分が物体領域を示し白い領域が空洞領域を示している。そして、仮想的な物理モデルにおいて、ｈ０＝０．１とし、ａ＝０．２とした。また、三角形要素を用いて離散化して有限要素法を用いて解析した。
　この２次元モデルにおいて得られた厚さ関数ｈ（ｘ）を図１０（Ｂ）に示す。
　また、特徴量算出において厚さ関数とは別の関数（スケルトン関数）を用いることにより、図１０（Ｃ）に示すように、形状のスケルトン（骨組）を算出することができる。
　その他、図１１（Ａ）に示す２次元の材料分布（２次元モデル）における幾何学的特徴量（厚さ、形状法線ベクトル、スケルトン）を、図１１（Ｂ）～（Ｄ）に示す。
　＜３次元モデル（板形状）の数値解析結果＞
　図１２（Ａ）に示す３次元の材料分布（３次元モデル）において幾何学的特徴量の算出の妥当性について検討した。
　この３次元モデルは、３次元空間において板形状をなすものであり、グレー部分が物体領域を示し、白い領域が空洞領域を示している。そして、仮想的な物理モデルにおいて、ｈ０＝０．１とし、ａ＝０．２とした。また、テトラ要素を用いて離散化して有限要素法を用いて解析した。
　この物理モデルにおいて得られた厚さ関数ｈ（ｘ）を図１２（Ｂ）に示しており、本実施形態により３次元の材料分布においても特徴量（厚さ）を算出できることが分かる。
＜２．本実施形態の効果＞
　このように構成した幾何学的特徴量算出装置１００によれば、線形偏微分方程式を用いた仮想的な物理モデルにより、画像データから構造物の幾何学的特徴量に関連する値を示す関連値データを算出し、当該関連値データから構造物の幾何学的特徴量を示す特徴量データを算出することができる。これにより、画像データにおける構造物の形状の境界を抽出すること無く、構造物の幾何学的特徴量を算出することができる。
＜３．本発明の変形実施形態＞
　なお、本発明は前記実施形態に限られるものではない。
　例えば、前記実施形態の幾何学的特徴量算出装置１００を用いて、医療用画像分析システムを構成しても良い。この場合、画像データ受付部１は、人の身体部位をスキャンして得られる医療用画像データ（例えばＣＴスキャンデータ）を受け付ける。そして、特徴量データ算出部３は、身体部位の幾何学的特徴量を示す特徴量データを算出する。
　また、前記実施形態の幾何学的特徴量算出装置１００を用いて、リバースエンジニアリング支援システムを構成しても良い。この場合、画像データ受付部１は、ハードウエア製品の構造を示す画像データを受け付ける。そして、特徴量データ算出部３は、ハードウエア製品の構造の幾何学的特徴量を示す特徴量データを算出する。
　さらに、前記実施形態の幾何学的特徴量算出装置１００を用いて、構造最適化システムを構成してもよい。この場合、画像データ受付部１は、構造最適化手法により演算装置が演算した構造物の画像データを受け付ける。ここで、構造最適化手法には、寸法最適化、形状最適化又はトポロジー最適化が考えられる。そして、特徴量データ算出部３は、構造物の幾何学的特徴量を示す特徴量データを算出する。
　この構造最適化システムであれば、最適化の過程において形状の寸法を定量的に知ることができるので、構造最適化手法の制約条件として用いることができる。その結果、最大形状寸法、最小形状寸法、又は形状寸法一定等の所定の寸法制約下で構造最適化を行うことができる。構造最適化への適用における最大の特長は、最適化の過程において、必ずしも制約を満たす必要がないので、初期設計解の依存性が低く、高い設計自由度を維持しながら、構造最適化が可能となる。
＜４．構造最適化システム＞
　次に、構造最適化システムの一例として、３Ｄプリンタ等に用いられる積層造形法により製造が可能な構造物を創成する構造最適化システムを説明する。
　この構造最適化システムは、レベルセット法を用いてトポロジー最適化を行うものである。なお、構造最適化システムは、前記実施形態の幾何学的特徴量算出装置１００が組み込まれたものであり、ＣＰＵに加えて揮発メモリやＨＤＤ等の記憶部を有する演算装置を備え、さらにマウスやキーボード等の入力手段、演算結果を出力するためのディスプレイやプリンタ等からなる出力手段を接続するための入出力インターフェイス等を有した汎用乃至は専用のコンピュータである
＜トポロジー最適化の考え方＞
　構造最適化問題とは、ある評価関数Ｊが最小（あるいは最大）となる構造物の形状Ω１を求める問題である。評価関数は、剛性、揚力、熱伝導特性などの物理的な特性を含む場合が多い。そのため、評価対象とする物理現象を表現する支配方程式を満足することを前提として、最適な形状を求める。支配方程式は、最適化問題における制約条件となり、基本的な構造最適化問題は次のように定式化される。
Figure JPOXMLDOC01-appb-M000015
　ここで、ｕは支配方程式の解として得られる状態変数、ｊ（ｕ）は所望の目的関数を表現する関数、Ｕａｄは、許容されるΩ１がなす集合を表す。評価関数Ｊに含まれる支配方程式の解ｕと形状Ω１の関係は、陽な形式の関数で表現されない場合が多く、支配方程式の数値解析解として与えられる。次に、構造最適化問題（数　）に対して、トポロジー最適化の適用について考える。まず、設計の対象とする構造物が配置しうる領域ΩＤ⊂Ｒｎを導入する。ただし、ｎは対象とする空間次元を表す。ここで、領域ΩＤは最適設計解の探索過程において変化しないため、固定設計領域と呼ばれる。固定設計領域は、構造により占められた領域Ω１（以下、構造領域）と構造に占められていない領域Ω０（以下、空洞領域）により構成されており、それらを識別する特性関数χを次式により定義する。
Figure JPOXMLDOC01-appb-M000016
　ただし、構造領域Ω１と空洞領域Ω０の境界は構造領域に含めるものとする。このとき、特性関数χを用いて、トポロジー最適化問題は次のように定式化される。
Figure JPOXMLDOC01-appb-M000017
　ただし，Ｆａｄは許容される形状の特性関数のなす集合を表す。トポロジー最適化問題では、設計変数が特性関数χであるから、最適構造設計問題が材料の最適分布問題として表現される特徴を持つ。その結果、設計変更の過程において、構造領域Ω１のトポロジーが変更される可能性を持つ。
＜トポロジー最適化問題の緩和と正則化＞
　トポロジー最適化問題は、最適解が存在しないため、不良設定問題である(Allaire, 2002)。従って、何らかの手段を用いて良設定問題に置き換える手続を必要とする。設計空間の緩和に基づく方法として、均質化設計法(Bendsoe and Kikuchi, 1988)が代表的である。また、広く普及している密度法(Bendsoe and Sigmund, 1999)は、均質化設計法の考え方に基づいて簡便な問題に置き換えた方法である。これらの方法では、構造領域と空洞領域の中間状態を許容するため、連続な密度関数を用いて最適構造が表現される。他方、正則化を用いた方法としては、レベルセット法に基づく方法(山田他,2009; Yamada et al., 2010)がある。この方法では，レベルセット関数と呼ばれるスカラー関数の等位面により構造領域の境界面を表現し、レベルセット関数の変更により構造領域の形状変更を表現する方法である。また、レベルセット関数に滑らかさを適切に確保することで、トポロジー最適化問題の正則化を行う。ここでは、外形形状の幾何学的制約について扱うため、レベルセット法に基づく方法を用いる。
＜レベルセット法に基づく方法＞
レベルセット法に基づく方法では、次式に示すスカラー関数φ（ｘ）を導入し、形状表現を行う。
Figure JPOXMLDOC01-appb-M000018
　レベルセット関数を用いて特性関数χを再定義すれば、次式となる。
Figure JPOXMLDOC01-appb-M000019
　レベルセット関数により定義された特性関数χを用いてトポロジー最適化問題を定式化する。
　次に、レベルセット関数の更新方法について述べる。基本的な考え方は。最適化の過程におけるレベルセット関数φ（ｘ）の変更を、仮想的な時間ｔにおける時間変化とすることである。レベルセット関数φ（ｘ）は、値域制約（－１≦φ（ｘ）≦１）を満たしながら、次式に示す反応拡散方程式に基づいて更新される。
Figure JPOXMLDOC01-appb-M000020
　ただし、Ｋ∈Ｒ＋は比例定数、Ｃ∈Ｒ＋は正規化定数、τ∈Ｒ＋は正則化定数、Ｊ’はトポロジー導関数(Amstutz et al.,2006)である。反応拡散方程式における拡散項は、レベルセット関数φ（ｘ）の滑らかさを保証するだけではなく、拡散項の相対的な寄与度を与える正則化係数τの調節により、得られる最適構造の幾何学的複雑さの定性的な制御ができる特徴も持つ。さらには、正則化係数の寄与を方向毎に異なる大きさに設定すれば，幾何学的複雑さに指向性を持たせることが可能となる。その結果、フライス加工等に要求される等断面形状制約を容易に考慮できる(Yamada et al., 2010)。また、比例定数Ｋと正規化定数Ｃは、レベルセット関数の値域制約が１程度であるから、同程度の大きさに設定するパラメータである。
＜幾何学的特徴量に対する制約を考慮したトポロジー最適化＞
＜幾何学的特徴量に対する偏微分方程式系＞
　幾何学的特徴量を抽出する対象形状Ω１を包含する領域Ωにおいて、次に示す偏微分方程式系を考える。
Figure JPOXMLDOC01-appb-M000021
　ここで，ｓｉ∈Ｈ１（Ω）はｘｉ方向に対する状態変数、ｅｉは標準基底、ａ∈Ｒ＋は拡散係数である。なお，拡散係数ａは相対的な影響が小さくなるように、無次元化拡散係数ａ～：＝ａ／Ｌ２を十分小さな値ａ～≪１に仮定する。ただし、Ｌは代表長さである。ここで、上記の偏微分方程式の力学的意味を議論するために、状態変数ｓｉから構成されるベクトル場ｓ＝［ｓ１ｓ２・・・ｓｎ］Ｔ及びその試験関数∀ｓ～∈Ｈ０ １（Ω）ｎを導入すると、次式の等価な弱形式が得られる。
Figure JPOXMLDOC01-appb-M000022
　ベクトル場ｓの支配方程式は、縦弾性係数がａ、ポアソン比が０である等方性線形弾性体で占められた領域Ωに対して、周囲∂Ωを完全拘束、領域Ω０に大きさが１の等方性ばねを配置、対象形状の境界∂Ω１から外向きに大きさ１の界面力を作用させた場合の弾性方程式である。従って、ベクトル場ｓは、このような力学的条件における変位場と等価である。直感的には、ポアソン比が０であるから、ベクトル場ｓの向きは。境界∂Ω１において、その法線ベクトルの向きと一致すると推察される。従って、対象形状Ω１の法線ベクトルを近似的に表現する関数ｆｎ（｛ｓｉ｝１≦ｉ≪ｎ；εｓｎ）は次式で与えられる。
Figure JPOXMLDOC01-appb-M000023
　ただし，εｓｎ∈Ｒ＋は、ｓ＝０の場合においてもベクトル関数ｆｎが定義できるように導入された微小な正の値を持つパラメータである。なお、ｓ＝０となるのは、境界面∂Ω１から十分に遠い場所であるから、微小パラメータεｓｎの設定値がベクトル関数ｆｎの意味に実質的に影響を与えないことを注記する。形状の境界面∂Ω１における厳密な法線ベクトルは、正則な境界面上のみで定義されるが、ベクトル関数ｆｎは、接ベクトルの勾配が不連続となる幾何学的特異点を含めて定義できる。従って、ベクトル関数ｆｎは、正則な境界面近傍においては、法線ベクトルに漸近する関数であり、本来の定義より広い意味で法線ベクトル定義し、領域Ω全域で値を持つ。証明には、偏微分方程式の解の振る舞いに関する厳密な議論が必要であり、今後の検討課題ではあるが、数値解析に基づく検証の結果、高い精度で近似されることが確認されている(Yamada, 2019b)。また、偏微分方程式における拡散項の特徴から、半径がＯ（√ａ）の範囲で状態変数ｓｉが平均化された振る舞いとなる。つまり、ベクトル関数ｆｎは、∂Ω１の近傍において、半径Ｏ（√ａ）の大きさで平均化された法線ベクトルであると言える。
　次に、状態変数ｓｉの分布に関する基本的な特徴を議論するために、２次元円形形状に対する偏微分方程式の解析解を示す。円の半径をＲとし、円の中心を原点とするとき、解析解は、極座標表示（ｒ，θ）を用いて、次式となる。
Figure JPOXMLDOC01-appb-M000024
　ただし、Ｋｍ（ｘ）は第二種変形ベッセル関数である。上式（数２３）に示すように、形状Ω１の内部においては、原点を通る線形関数、外部では急激かつ単調に減小しながら零に収束する関数である。したがって、ベクトル場ｓの大きさの絶対値は。境界∂Ω１に近づくにつれて単調に増加する。また、ベクトル場ｓの向きは、境界∂Ω１において、境界の法線方向と厳密に一致している。領域Ω１の内部と外部においても、最も近い境界面の法線方向と一致している。
＜オーバーハング制約に対する幾何学的制約条件＞
　次に、積層造形法におけるオーバーハング形状に関する幾何学的制約条件を明確化する。最初に、図１３（ａ）に示すように、造形方向と造形対象の形状（以下、造形形状）の接ベクトルとのなす角を積層傾斜角θａと定義する。積層傾斜角θａが一定値を超えると、造形中の自重の作用により、造形物の表面精度が著しく悪化する問題がある。例えば、図１３（ａ）の場合、破線で示した造形形状は、積層造形法による製造が困難であると言える。このように、下部に十分な支持がないために製造が困難な部分形状をオーバーハング形状、その部分形状の水平面からの傾斜角θをオーバーハング角と呼ぶ。この問題を回避するために、オーバーハング形状の下部に、サポート材と呼ばれる仮の造形物を配置する対処法が用いられる。しかしながら、サポート材の配置に基づく対処法は、製造工程やコストの増加を招くだけではなく、サポート材の除去が可能な単純な造形形状に限定される。
　そこで、サポート材なしに積層造形を可能とする造形形状の幾何学的条件を明確化する。先行研究の多くは、形状の法線ベクトルに基づいて幾何学的制約条件を定式化している。法線ベクトルに基づく考え方を図１３（ｂ）に示す造形形状に適用すれば、破線で示した境界が制約違反と判断される。しかしながら、点Ｐｓ １＋と点Ｐｓ ２＋及びその上部は、明らかに積層造形が不可能な形状である。
　本実施形態では、造形が成立しない領域の評価及び平均曲率に基づいた幾何学的特異点の評価に基づく方法を提案する。具体的には、図１４（ａ）に示す円形形状において，積層方向をｘ２軸方向と設定した場合、境界ｌｃを評価するのではなく、境界ｌｃが円弧となる扇形の領域は、制約違反により成立する領域（以下、制約違反領域）として定義する。この領域は、図１４（ｂ）に示す左側への積層造形ができない領域（以下、負側違反領域）と、図１４（ｃ）に示す右側への積層造形ができない領域（以下、正側違反領域）のいずれにも該当する領域である。なお、３次元の場合は、法線ベクトルを水平面へ投影した場合の向き及び鉛直方向が基底ベクトルとなる２次元平面を考え、同様の議論を考えれば良い。
　次に、図１３（ｂ）の点Ｐｓ １＋と点Ｐｓ ２＋のような幾何学的特異点の抽出法について考える。ここで、点Ｐｓ １＋と点Ｐｓ ２＋のような凸型形状の幾何学的特異点では、平均曲率は正の無限大であり，点Ｐｓ １－と点Ｐｓ ２－のような凹型形状の幾何学的特異点では平均曲率は負の無限大であることに注目する。平均曲率は法線ベクトルの発散により与えられるから、幾何学的特異点の近傍領域は、ベクトル関数ｆｎの発散の絶対値が大きい領域であり、その正負により凸型形状もしくは凹型形状であるかを識別できる。なお、前述のように、ベクトル関数ｆｎはＯ（√ａ）程度の範囲で平均化されるため、幾何学的特異点近傍においても、ベクトル関数ｆｎの発散は有界な値を持つことを注記する。
　以上の議論から、積層造形におけるオーバーハング形状に対する幾何学的制約条件は、制約違反領域及び制約違反領域近傍の凸型形状の幾何学的特異点が空集合となるように定式化すれば良い。
＜幾何学的制約条件の定式化＞
　幾何学的特徴量に対する偏微分方程式系及び上記の議論に基づき、２次元問題における幾何学的制約条件を定式化する。単純化のために、積層方向がｘ２の正の方向となるように座標系を選ぶこととする。このとき、積層限界方向を与える単位ベクトルξ１及びξ２は、幾何学的関係から次のように与えられる。
Figure JPOXMLDOC01-appb-M000025
 　負側違反領域は、ベクトルξ１と拡張した法線ベクトルｆｎの内積が正となる対象形状Ω１の内部であるから、ベクトル場ｓは対象形状Ω１近傍のみで非零であることに注目すれば、領域Ω全体における負側違反領域近傍は、ベクトルξ１とベクトル場ｓの内積が正の領域として与えられる。同様に、正側違反領域は、ベクトルξ２とベクトル場ｓの内積が正となる領域として与えられる。従って、制約違反領域に基づいた幾何学的制約条件は次式により与えられる。
Figure JPOXMLDOC01-appb-M000026
　ただし、Ｒはランプ関数である。なお、ランプ関数の代わりにヘビサイド関数を用いた定式化も可能であるが、ヘビサイド関数の導関数はデルタ関数であるから、最適化問題の多峰性を強めてしまう問題を持つことを注記する。さらに、ベクトル場ｓはＯ（１／√ａ）の大きさを持つことに注目して適切に正規化を行うと、幾何学的制約条件は次のよう定式化される。
Figure JPOXMLDOC01-appb-M000027
　ここで，Ｔ（ｘ；γ）は次式である。
Figure JPOXMLDOC01-appb-M000028
　パラメータγを十分に大きく設定した場合、上式の被積分関数値は、制約違反領域の内部において一定（Ｔ＝１）となる。なお、空間分布を正規化するＴは、次に示す四つの条件を満たす必要がある。
Figure JPOXMLDOC01-appb-M000029
　次に、幾何学的特異点について考える。前述のように、拡張された法線ベクトルの発散の絶対値が大きい制約違反領域近傍が幾何学的特異点の近傍である。さらに、その符号が正の領域は凸型形状であるから、制約汎関数値Ｇ＾ａｍの被積分関数を過大に評価すれば良い。同様に符号が負の領域は凹型形状であるから、制約汎関数値Ｇ＾ａｍの被積分関数を過小評価すれば良い。ただし、負の値をとることを許容してはならない。従って、幾何学特異点を考慮した制約汎関数は次式となる。
Figure JPOXMLDOC01-appb-M000030
 　ここで、β∈Ｒ＋は、幾何学的特異点の排除度合いを与えるパラメータである。パラメータβの設定値が大きい場合、制約違反に近い形状、かつ、局所的に曲率が大きい場所を排除する。そのため，許容されるオーバーハング角の制約を違反する部分構造を許容する。従って、幾何学的特異点に対する制御パラメータβは、下向き凸型形状の幾何学的特異点を排除しながら、適切な曲率を維持できる範囲に設定する必要がある。パラメータβの設定値の具体的な影響については、後述の数値解析例で示す。
＜最適化問題の定式化＞
　次に、弾性テンソルがＤの等方性線形弾性体で構成される構造領域Ω１について考える。構造領域Ω１は、全境界∂Ωの一部の境界Γｕで変位拘束、別の一部の境界Γｔで表面力ｔ∈Ｌ２（Γｔ）２が負荷されているとする。このとき、変位場をｕ∈ＨΓｕ １（Ω１）２とし、平面ひずみを仮定された二次元問題における体積制約及びオーバーハング制約を考慮した平均コンプライアンス最小化問題は次のように定式化される。
Figure JPOXMLDOC01-appb-M000031
　ただし、ｈは領域の厚み、Ｇｍａｘ∈Ｒ＋は許容される体積の上限値、ε（ｕ）：＝１／２×（∇ｕ＋∇ｕＴ）はひずみであり、ＨΓｕ １は次に示す関数空間である。
Figure JPOXMLDOC01-appb-M000032
　また，弾性テンソルＤは、縦弾性係数Ｅとポアソン比ｖを用いて次式で与えられる。
Figure JPOXMLDOC01-appb-M000033
　ただし、δｉｊはクロネッカーのデルタである。
＜数値実装法＞
＜最適化アルゴリズム＞
　構造最適化システムの演算装置は、次に示すアルゴリズムに従って最適化を行う。
（ステップ１）
　レベルセット関数の適当な初期値φ０を設定する。
（ステップ２）
　有限要素法を用いて、変位場及びベクトル場ｓを数値解析する。
（ステップ３）
　目的汎関数及び制約汎関数を計算し、収束していれば最適化を終了する。
（ステップ４）
　有限要素法を用いてベクトル場ｓに対する随伴場を数値解析し、随伴変数法により設計感度（トポロジー導関数）を求める。
（ステップ５）
　反応拡散方程式を用いてレベルセット関数φを更新し、（ステップ２）へ戻る。
　以上の手順により最適構造の候補を示すレベルセット関数値を得る。
　上記（ステップ１）のレベルセット関数は、演算装置のレベルセット関数データ格納部に格納されている。レベルセット関数は、初期構造が設定された設計領域の各部が、構造を形成し、物体により占められた物体領域（物体相）、空洞を形成する空洞領域（空洞相）、又はそれら領域の境界であるかを示すものであり、物体領域を表す値及び空洞領域を表す値の間の所定値が、物体領域及び空洞領域の境界を表す。なお、レベルセット関数データは、例えばユーザが入力手段を用いることにより入力される。
　また、設計領域を示す設計領域データは、演算装置の設計領域データ格納部に格納されている。なお、設計領域（当該設計領域を要素分割する構造格子（メッシュ）情報を含む。）を示す設計領域データは、例えばユーザが入力手段を用いることにより入力される。
　さらに、前記設計領域（固定設計領域）の境界条件を示す境界条件データは、演算装置の境界条件データ格納部に格納されている。具体的な境界条件としては、例えば設計領域の拘束条件、初期構造に作用する荷重等の外力（表面力）などである。なお、境界条件データは、例えばユーザが入力手段を用いることにより入力される。
　その他、最適構造を求めるときに用いる体積制約に必要な制約の条件値及び変位場の解析に必要なヤング率やポアソン比の値等の材料定数である解析データは、演算装置の解析データ格納部に格納されている。なお、解析データは、例えばユーザが入力手段を用いることにより入力される。
　そして、（ステップ２）～（ステップ５）は、演算装置のレベルセット関数更新部により実行される。レベルセット関数更新部は、所定の制約条件下において、剛性、固有振動数等の構造物の性能を目標値に近づけるようにレベルセット関数を更新して、物体領域と空洞領域との境界を移動させるとともに、レベルセット関数の更新に伴う物体領域内のトポロジー変化（形態変化）を許容して物体領域内に新たな空洞領域（穴）を形成し、その新たな空洞領域と物体領域との境界を移動させるものである。
＜オイラー座標系に基づく変位場の近似解法＞
　構造領域Ω１はオイラー座標系により表現されており、最適化の過程における反復計算ごとに、有限要素の生成が必要となる。計算コストの観点から、Ersatz material approach（Allaire et al., 2004）による近似解法を適用する。すなわち、空洞領域Ω０は相対的に小さい弾性テンソルを持つ構造材料と仮定し、その界面近傍では滑らかに分布する材料特性を持つと仮定する。具体的には、有限要素解析において、次に示す拡張した弾性テンソルＤ～を用いて、変位場に対する支配方程式を固定設計領域ΩＤに対して拡張する。
Figure JPOXMLDOC01-appb-M000034
　ただし、εｕ∈Ｒ＋は、近似ヘビサイド関数Ｈ～の遷移幅を与える微小パラメータ、ｄ∈Ｒ＋は、構造領域の縦弾性係数に対する空洞領域の相対値である。
＜ベクトル場ｓに対する特性関数の拡張＞
　ベクトル場ｓは固定設計領域ΩＤを包含する領域Ωにおいて解析する必要があるため、領域Ωへ広げた拡張平滑化特性関数Ψ∈Ｈ１（Ω）ｎを次式により定義する。
Figure JPOXMLDOC01-appb-M000035
　ただし、εΨ∈Ｒ＋は、平滑化度合いを与える微小パラメータ、εｗ∈Ｒ＋は、特性関数の遷移幅を与える微小パラメータである。
　ランプ関数Ｒ（ｘ）はｘ＝０において強い意味で微分が定義できないため、微小パラメータεｓｒ∈Ｒ＋を導入し、ランプ関数Ｒに漸近する双曲線を近似ランプ関数Ｒ～として定義する。
Figure JPOXMLDOC01-appb-M000036
 　設計感度を導出する感度解析の段階において、近似ランプ関数を用いて数値解析を行う。
＜数値解析例＞
　数値解析に用いたモデルは、図１５（ａ）に示すように、６×１の矩形領域を固定設計領域とし、この下辺と同一直線上に積層プラットフォームを設定し、積層方向は鉛直上向きとする。このとき、縦弾性係数が１、ポアソン比が０．３の平面ひずみを仮定された等方性線形弾性体で占められた形状Ω１の体積制約付き剛性最大化問題について考える。
　力学的条件は、図１５（ｂ）に示すように、下辺両端の境界Γｕにおいてローラー支持し、上辺中央下向きに大きさｔ＝１の表面力を作用させる、左右対称であるから、中央に対称条件を与えることで、右半分の領域のみを解析する。ベクトル場ｓを定義する領域Ωは、固定設計領域を包含するように、左辺、右辺及び上端側に幅０．２を持つ矩形領域分を広げた領域とした。下辺側に超過した制約違反領域を評価する必要がないため、領域を広げる必要はないことを注記する。また、許容される体積の上限値は固定設計領域の５０％する。固定設計領域は、１５０×５０の構造格子に対して、格子中心周りに対称４分割した３角形要素を用いて有限要素分割し、変位場、幾何学状態場及びレベルセット関数場を１次有限要素により離散化する。反応拡散方程式及び幾何学的特徴量に対する偏微分方程式の無次元化に用いる代表長さＬはＬ＝１とする。また、無次元化拡散係数ａ～は幾何学的特徴量の評価をする際の平滑化半径の大きさを与えるから、幾何学的複雑さを与える正則化係数τと比較して十分に小さく設定する必要がある。ここでは、ａ～＝０．０３τと設定する．オーバーハング制約汎関数Ｇａｍの正規化パラメータγはγ＝２．５、ランプ関数の平滑化パラメータεｓｒは、適切に滑らかさを保持できるように、ξ・ｓに対してはεｓｒ＝３、曲率評価（１＋β×ｄｉｖｆｎ）に対してはεｓｒ＝１×１０－９とする。変位場の近似解法に用いるパラメータは、ｄ＝１×１０－３、εｕ＝０．５とする．幾何学的特徴量に対する偏微分方程式の近似解法に用いるパラメータは、εΨ＝１×１０－４Ｌ２、εｗ＝０．９５とする。
　最初に、オーバーハング制約を考慮しなかった場合の最適構造を図１６に示す。図中、黒色で占められた領域が構造領域Ω１、灰色で絞められた領域が空洞領域Ω０を表す。図１６に示すように、幾何学的正脚条件を考慮しない場合、積層フラットフォームと平行に近い部分構造を持つ最適構造が得られる。また、正則化係数τは構造の幾何学的複雑さを制御する効果を持つものの、正則化係数τの調整を行っても、積層造形が困難な部分構造を排除できないことを確認できる。
＜幾何学的特異点に対する制御パラメータβの影響＞
　幾何学的特異点に対する制御パラメータβが最適構造に与える影響について検討する。許容される最小のオーバーハング角をθ＝４５°に設定し、幾何学特異点に対する制御パラメータβをβ＝３×１０－２、２×１０－２、１×１０－２、８×１０－３、３×１０－３、５×１０－４、１×１０－４、０に設定した。得られた最適構造を図１７に示す。β＝０の場合、幾何学的特異点に対する考慮がないため、図１７（ｈ）に示すように、凸型の幾何学的特異点を持つ部分構造が得られている。制御パラメータβの設定値が大きい場合、図１７（ａ）に示すように，許容されるオーバーハング角の制約を違反する部分構造の領域が広いことが確認できる。制御パラメータβが小さくなるにつれて、許容される曲率の範囲が広くなるため、凹型の幾何学的特異点に近づいていくことも確認できる。従って、制御パラメータβを適切な範囲に設定すれば、オーバーハング制約を満たさない部分領域が非常に小さく、積層造形が可能な最適構造が得られることがわかった。また、幾何学的特異点を考慮しなければ、所望のオーバーハング制約を満たす最適構造が得られないこともわかった。
　なお、以降の数値解析例においては、β＝６×１０－３に設定する。
＜正則化係数τの影響＞
　次に、異なる正則化係数τを設定し、最適化を図った。許容される最小のオーバーハング角はθ＝４５°とする。このとき、得られた最適構造を図１８に示す。図１８に示すように、いずれの構造も、所望のオーバーハング制約を満たす構造が得られていることが確認できる。また、正則化係数τの設定値の大きさにより、最適構造の幾何学的な複雑さが異なることを確認できる。従って、オーバーハング制約を考慮した場合においても、最適構造の幾何学的複雑さを定性的に制御可能であることがわかった。
＜最小許容オーバーハング角の影響＞
　最後に、許容される最小オーバーハング角の設定値が最適構造に与える影響について検討する。許容される最小のオーバーハング角として、４つの設定値（θ１５°、３０°、４５°、６０°）について最適化を図った。正則化係数は、τ＝３×１０－３及びτ＝４×１０－４に設定し、得られた最適構造を図１９に示す。図１９に示すように、いずれの結果も、オーバーハング制約を考慮した結果が得られ、オーバーハング角の設定値毎に異なる最適構造が得られた。
　上記の構造最適化システムは、オーバーハング形状に関する幾何学的制約条件を加味したトポロジー最適化を行うものであったが、幾何学的制約条件としては、オーバーハング形状に限られない。例えば、図２０に示すモデルを用いた体積制約付き剛性最大化問題において、構造物の厚さ（幅）を所定値以下とする幾何学的制約条件を加味したトポロジー最適化を行うようにしてもよい。その結果を図２１に示す。
　その他、本発明は前記実施形態に限られず、その趣旨を逸脱しない範囲で種々の変形が可能であるのは言うまでもない。
１００・・・構造物の幾何学的特徴量算出装置
１　　・・・画像データ受付部
２　　・・・関連値データ算出部
３　　・・・特徴量データ算出部
４　　・・・特徴量データ出力部

Claims (10)

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

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2020518248A JP7224675B2 (ja) 2018-05-07 2019-04-24 構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、幾何学的特徴量算出プログラム

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2018089391 2018-05-07
JP2018-089391 2018-05-07

Publications (1)

Publication Number Publication Date
WO2019216221A1 true WO2019216221A1 (ja) 2019-11-14

Family

ID=68468016

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2019/017419 WO2019216221A1 (ja) 2018-05-07 2019-04-24 構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、幾何学的特徴量算出プログラム

Country Status (2)

Country Link
JP (1) JP7224675B2 (ja)
WO (1) WO2019216221A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2022105478A (ja) * 2021-01-04 2022-07-14 国立大学法人 東京大学 構造物の設計支援装置、処理対象物の処理支援装置、設計支援方法および設計支援プログラム

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002007487A (ja) * 2000-06-27 2002-01-11 Asahi Kasei Corp 最適形状の設計方法及びこれを用いた最適形状の設計システム
WO2008023491A1 (fr) * 2006-08-24 2008-02-28 Olympus Medical Systems Corp. dispositif et procédé de traitement d'images à usage médical
WO2010029810A1 (ja) * 2008-09-11 2010-03-18 国立大学法人京都大学 構造最適化装置、構造最適化方法及び構造最適化プログラム

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2478329B (en) 2010-03-03 2015-03-04 Samsung Electronics Co Ltd Medical image processing

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002007487A (ja) * 2000-06-27 2002-01-11 Asahi Kasei Corp 最適形状の設計方法及びこれを用いた最適形状の設計システム
WO2008023491A1 (fr) * 2006-08-24 2008-02-28 Olympus Medical Systems Corp. dispositif et procédé de traitement d'images à usage médical
WO2010029810A1 (ja) * 2008-09-11 2010-03-18 国立大学法人京都大学 構造最適化装置、構造最適化方法及び構造最適化プログラム

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2022105478A (ja) * 2021-01-04 2022-07-14 国立大学法人 東京大学 構造物の設計支援装置、処理対象物の処理支援装置、設計支援方法および設計支援プログラム
JP7365718B2 (ja) 2021-01-04 2023-10-20 国立大学法人 東京大学 構造物の設計支援装置、処理対象物の処理支援装置、設計支援方法および設計支援プログラム

Also Published As

Publication number Publication date
JP7224675B2 (ja) 2023-02-20
JPWO2019216221A1 (ja) 2021-05-13

Similar Documents

Publication Publication Date Title
Jahangiry et al. An isogeometrical approach to structural level set topology optimization
Kang et al. Isogeometric topology optimization of shell structures using trimmed NURBS surfaces
Nguyen et al. Three-dimensional topology optimization of auxetic metamaterial using isogeometric analysis and model order reduction
Philipp et al. Integrated design and analysis of structural membranes using the Isogeometric B-Rep Analysis
Dilgen et al. Topology optimization of acoustic mechanical interaction problems: a comparative review
Wang et al. Structural shape and topology optimization in a level-set-based framework of region representation
Ye et al. Topology optimization of conformal structures on manifolds using extended level set methods (X-LSM) and conformal geometry theory
Seo et al. Shape optimization and its extension to topological design based on isogeometric analysis
Kruse et al. Isogeometric collocation for large deformation elasticity and frictional contact problems
EP3092624B1 (en) Method for structure preserving topology optimization of lattice structures for additive manufacturing
Wang et al. A level set method for structural topology optimization
Lawry et al. Level set shape and topology optimization of finite strain bilateral contact problems
Da Silva et al. Stress-based topology optimization of continuum structures under uncertainties
Hildebrandt et al. Interactive surface modeling using modal analysis
Kang et al. Isogeometric shape optimization of trimmed shell structures
Liu et al. Topology optimization of structures considering local material uncertainties in additive manufacturing
Stutz et al. Singularity aware de-homogenization for high-resolution topology optimized structures
Kedward et al. Gradient-limiting shape control for efficient aerodynamic optimization
US11501032B2 (en) Surface developability constraint for density-based topology optimization
Kumar et al. Synthesis of c0 path-generating contact-aided compliant mechanisms using the material mask overlay method
Ullah et al. A boundary element and level set based topology optimisation using sensitivity analysis
Tur et al. Stabilized method of imposing Dirichlet boundary conditions using a recovered stress field
WO2019216221A1 (ja) 構造物の幾何学的特徴量算出装置、医療用画像分析システム、リバースエンジニアリング支援システム、構造最適化システム、及び、幾何学的特徴量算出プログラム
Kerfriden et al. Statistical extraction of process zones and representative subspaces in fracture of random composites
Riehl et al. On the discrete variant of the traction method in parameter-free shape optimization

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19798965

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2020518248

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19798965

Country of ref document: EP

Kind code of ref document: A1