Paper
An Optimal Model Subdivision Method for 3D FDM
Printers
Satoshi Nakajima † , Mikio Shinya (member)††,††† , Takashi Yuasa †† ,Michio Shiraishi ††Abstract
Although FDM 3D printers have been getting popular, it is still difficult for novice users to determine basic parameters suchas deposition directions and model subdivision. This paper describes an optimization method that determines the optimalmodel subdivision and deposition direction. The cost function is based on the supporting structure volume and the numberof regions to adhere, which can be computed from multi-layered depth maps on GPU. The optimization can be efficiently
achieved by the dynamic programing. Experiments showed that the method can improve object surface smoothness andreduce working time of support removal.

Key words: 3D printer, optimization, FDM, dynamic programing.

1. Introduction
A 3D printer is a device that creates 3D physical objects from geometric models. In recent years,
low-cost consumer-level 3D printers have been released
and rapidly spread out. Many attractive applications
have been developed, such as animatronic mold1)2) andKinect-base self-portrait printing3) . This results inmore opportunities for non-professional novice users tocreate 3D objects by 3D printers, and it is important todevelop tools that assist untrained users to easily createhigh quality objects.
Zhou developed a structural analysis method that
identifies structural mechanical problems4) . Reductionof used materials was discussed by Wang5) , who proposed to construct inner objects by truss scaffoldings
and to apply nonlinear optimization. Prévost et. al.,proposed a method balancing static models at rest6) .
However, many issues still remain open.
Most of low-cost systems adopt the fused deposition
modeling (FDM) method, which makes physical objects by depositing melted materials such as ABS resins.
FDM 3D printers build up 3D objects layer by layer.
Therefore, the deposition direction has a large influenceThis paper reported in IWAIT2014, January 6 2014.
Received Month xx, 20xx; Revised Month xx, 20xx; Final received
Month xx, 20xx; Accepted Month xx, 20xx
† Fujitsu Broad Solution & Consulting, Inc.
(Tokyo, Japan)†† Graduate School of Science, Toho University(Chiba, Japan)††† UEI Research(Tokyo, Japan)

on the production. When FDM 3D printers make objects containing holes or bridge-structures, they have todeposit supporting structures underneath the objects,
as shown in Fig. 1-a. It is necessary to remove these
supports by hand after deposition, which often damages the objects, as shown in Fig. 1-b. Therefore, it isvery important to control support structures to obtainwell-finished objects and to reduce work time.
The support volume can be reduced by changing thesize or structure of supports7)8) . However, it is possible
to reduce support volumes by selecting an appropriate
deposition direction, as well. We can also reduce it bydividing a geometric model into several parts, creating each part separately, and adhering them after thedeposition, as shown in Fig. 2. For novice users, however, it is difficult to determine the appropriate direction and model subdivision. Therefore, it is important
to develop a computation model that calculates optimal
deposition directions and model subdivision.
Although a model subdivision technique based onconcave edges has been reported9) , it is restricted tosimple geometric objects consisting of polygons parallel to each other and it is difficult to apply the method
to complicated shapes such as figure dolls. To the authors’ knowledge, optimization methods to determine
model subdivision based on the cost function has notbeen discussed.
This paper describes an optimization method that determines the optimal model subdivision and deposition
direction by utilizing CG techniques and the dynamic
programing. The cost function we used in this paper

ITE Transactions on Media Technology and Applications Vol. xx, No. xx, pp. 1–7 (20xx)(1) 1was defined by a sum of supporting structure volumes
and the number of regions to adhere, although it canbe easily introduce other factors such as the area ofthe regions and weighting functions of support structure length.
Supporting structure volumes can be efficiently computed from multi-layered depth maps on GPU by using the depth peeling method10) . The number of adhesion regions is also calculated from the layered depth
maps using a filling algorithm. The optimization is also
efficiently achieved by the dynamic programing. We
restricted model subdivision planes to be perpendicular to the deposit direction, which allows us to apply
the dynamic programing. By using the optimal subdivided object, even novice users can reduce the problems
caused by the support removal.

(a)Supporting structures

Fig. 1

An example of supporting structures.

(a) Model subdivision.

Fig. 2(b) After their removal.

(b) After the adhesion.

Model subdivision and adhesion.

2. 3D Printer2. 1 3D printing systems
A 3D printer is a device creating 3D physical objects from geometric models. There are several types of
3D printing systems, such as the stereo-lithography, theselective laser sintering (SLS) and the fused deposition
modeling (FDM). Stereo-lithography 3D printers scana high power laser beam and draw the cross-sectionsof 3D physical objects onto the surface of the pool oflight curing resin, repeating it layer by layer. SLS 3D
printers create 3D physical object by fusing small particles such as nylon using a high power laser. FDM
3D printers create objects by depositing thermoplas2 (2)tic resin such as ABS resins or polycarbonates, by depositing the melted materials layer by layer. The standon which objects are placed horizontally moves to construct object cross-sections like a plotter and verticallysteps down after the completion of a layer. The FDM
system has enabled low-cost products and is most popular in the consumer-level market.
2. 2 Supporting Structures
When FDM 3D printers create objects with holesor bridge-structures, they have to deposit supporting
structures underneath the objects (Fig. 1-a). After thedeposition, supporting structures have to be removed
by hand, which often damages the objects. Therefore,
it is very important to control support structures toobtain well-finished objects and to reduce work time.
2. 3 Cost function
To determine appropriate directions and model subdivision, we have to consider the following factors.
( 1 ) Object qualities: The object is often damaged
by removing supporting structures. It can be expected that less support structures decrease thedamage.
( 2 ) Support removal time: Support removal byhand is a labor intensive task. This factor also
depends on the amount of supporting structures.
( 3 ) Adhering time: It is necessary to adhere divided parts, which also requires manual operations. This increases with respect to the numberof subdivision planes.
The volume of supporting structures and the adhering region can be considered to have dominant influences on the factors (1), (2) and (3). Thus, we adopted
a cost function consisting of two terms: the volume of
supporting structures and the number of regions to adhere after the deposition. Note that the volume termmonotonously decreases with the number of subdivision
while the number of regions to adhere increases.
The area of adhering regions can be another important factor. It will be also shown that this factor canbe also integrated to our cost function. Design issues
on the cost function will be discussed in Section 5.4.

3. Computation of Supporting Structure
Volumes
A supporting structure volume can be regarded asa volume of shadow regions when light rays come inthe deposit direction and can be calculated by using
depth values as Fig. 3. Let us assume model subdivision planes to be perpendicular to the deposit direction

ITE Transactions on Media Technology and Applications Vol. xx, No. xx (20xx)

maps. Referring to these maps, the supporting structure volume can be efficiently calculated.
In 3D printing, it is also possible to depose objects
in the reverse direction as well. Therefore, we set the
supporting structures volume Hs is the smaller volumes
in the direction s and −s, as:Z
Depositing direction0
Hs (D0 , D1 ) = min(Hs0 , H−s).

Fig. 3

The area of adhering regions, Ha , can be also calculated in a similar way. By setting Li = 1, the area Ha
is given by counting the number of points as:
X
X
1
Ha =
(7)
X outside object

Supporting structures.

Z
D1

(6)

The number of regions to adhere, M (z), can be calculated from the cross-sections of the objects and subdivision planes. A bit map image of the cross-section
CD (x, y) can be constructed from the multi-layered
depth map. The number of regions is then counted
through simple region filling.

z3
z2
z1

D0

4. Optimization Method
The cost function gs is defined by the summation of
the supporting structure volume Hs (Zi , Zi+1 ) and thenumber of regions to adhere M (Zi ) as:

X
Fig. 4

Length of supporting structure at X.

z. In Fig. 4, values zi (i = 1, . . . , n) indicate the depth
value of the object surface at X, and D1 and D0 are
the depth values of the upper subdivision plane and thelower plane:
D0 < zl < zl+1 < · · · < zh < D1 .
In this example, we have to set a support between z1
and z2 at X and its length is z2 − z1 . By summingup these support length, the total support length at X,
hs (X), can be calculated by:hs (X) =

X

Li ,

(1)

outside object
L0 = zl − D0 ,(2)

Li = zl+i+1 − zl+i ,(3)

Lh = 0.

(4)Note that Lh is zero because no support is necessary
above the top surface. The support volume in this subdivision section can be calculated by the sum of hs (X):X
Hs0 (D0 , D1 ) =hs (X).
(5)
X

Depth values zi (X) can be stored in multi-layered depth
Paper

gs (Zi , Zi+1 ) = Hs (Zi , Zi+1 ) + αM (Zi ),(8)

where α denotes a weight constant balancing the twoterms.
Selecting the optimal α value is generally a difficult
task. However, it can be roughly estimated by regarding it as a balancing factor between the dynamic rangesof a decrease in the support volume Hs and the numberof regions to adhere N . Given the expected decrease in
the volume δH at the maximum allowable number ofthe planes Nm , α can be estimated by δH/Nm . In thecase of Fig. 9, δH is about 0.15. If we set Nm = 15,
α = 0.15/15 = 0.01 would be a good value.
4. 1 Optimization of subdivision planes
Assume that an object is subdivided by N planes
perpendicular to the deposition direction s, as shownin Fig. 5. Let Zi (i = 1, · · · , N ) be their depth values,
and the minimum and maximum depth value of the object be Zo and Ze = ZN +1 . We quantize z-value innz levels. The total cost function, G, is defined by thesum of the cost values between neighboring subdivision
planes, gs (Zi , Zi+1 ), as:G=NX

gs (Zi , Zi+1 ).

(9)

i=0

(3) 3

Ze

s

㻰㼑㼜㼛㼟㼕㼠㻌㼐㼕㼞㼑㼏㼠㼕㼛㼚㻌㼟㼖㻌H㻯㼞㼑㼍㼠㼑㻌㼙㼡㼘㼠㼕㼘㼍㼥㼑㼞㼑㼐㻌㻌㼐㼑㼜㼠㼔㻌㼙㼍㼜㻌

Z3

㻯㼍㼘㼏㼡㼘㼍㼠㼑㻌㼠㼔㼑㻌㻴㼟㻔㼆㻘㼆䇻㻕㼍㼚㼐㻌㻹㻔㼆㻕㻌㻌

Z2
㻯㼍㼘㼏㼡㼘㼍㼠㼑㻌㻰㻼㻙㼠㼍㼎㼘㼑㻌

Z1

㻯㼍㼘㼏㼡㼘㼍㼠㼑㻌㻳㼙㼕㼚㻔㻺㻧㼟㼖㻕㻌㼒㼛㼞㻝䍺㻺䍺㻺㼙㼍㼤㻌H

Z0㻻㼜㼠㼕㼙㼕㼦㼑㼐㻌㼐㼑㼜㼛㼟㼕㼠㻌㼐㼕㼞㼑㼏㼠㼕㼛㼚㻌㼍㼚㼐㻌㼟㼡㼎㼐㼕㼢㼕㼟㼕㼛㼚㻌㼜㼘㼍㼚㼑㼟㻌
Fig. 5

Supporting structure volume.

Fig. 6The minimization of G can be efficiently achievedby the dynamic programing at an O(nz 2 ) computation
cost, while the cost for the exhaustive search is O(nNz ).
First, we calculate gs (z, z 0 ) and store the values into atable. Next, we calculate Ej (z) defined by

E0 (Z1 ) = gs (Z0 , Z1 ),Ej (Zj ) =(10)

min (Ej−1 (Zj−1 ) + gs (Zj−1 , Zj )).

Zj−1 <Zj

and also stored as DP-tables. The minimum values
Gmin for given N and S can be found by a linear search,as:

Gmin (N ; s) =min G

Z1 ,···,ZN= min(EN −1 (ZN ) + gs (ZN , Ze )). (11)ZN

Outline of the proposed method

We also construct a bitmap image of cross-sections
Cz (x, y) and count the number of adhesion, M (Z).
Using Hs (Z, Z 0 ) and M (Z), the DP-table Ek is calculated by Eq. (11). Minimal cost values Gmin (N ; sj )are then calculated form Eq. (11) for each number ofsubdivision planes N . We adopt N that provides thesmallest cost value. We execute this process for all sample directions and the optimal deposition direction is
the sample direction in which the cost value is smallest.
The object is clipped by the optimal subdivision
planes and the clipped polygons are triangulated. We
used Delaunay triangulation and implemented it using
CGAL11) .
In the experiments, we used multi-layered depth
maps at 64 × 64 resolution. We quantized depth valuesby 64 levels. We sampled 19 deposition directions asshown in Fig. 7.

4. 2 Optimization of number of subdivision
The optimization represented by Eq. (11) is conducted for each number of subdivision plane N =1, 2 · · · , Nmax . We adopt N that provides the smallest
cost value as the optimal number of planes.
4. 3 Optimizing of deposit direction
The cost function Eq. (9) is non-linear with respect
to the deposit direction. We decided to simply sample
deposition directions and calculate the minimal cost ineach direction. We adopt the direction that provides
the smallest cost value as the optimal deposit direction.

5. Experiments
5. 1 Procedure and implementation
Fig. 6 outlines the proposed method. We sample
deposit direction sj and calculate multi-layered depth
maps zi [x][y] by depth peeling method10) . By referringto the maps, we calculated support volumes Hs (Z, Z 0 ).
4 (4)Fig. 7

Sampling direction.

5. 2 Optimization
We applied the proposed method to the Stanford

ITE Transactions on Media Technology and Applications Vol. xx, No. xx (20xx)

Buddha model12) and made physical objects. In thisexperiment, we fixed the deposition direction and thenumber of subdivision N = 1 to show the basic properties.
Fig. 8 shows the cost value G when we changed theweighting factor α in Eq. (8). The supporting structure
volume is normalized by the volume of the object. As
shown in the figure, G has the minimum value near thecenter of the object.

㻝㻚㻞
1.5

㻝

#16
1

㻳㻌

㻜㻚㻤

ߙ ൌ ͲǤͲʹ

㻜㻚㻢
㻜㻚㻠

0.5

ߙ ൌ ͲǤͲͳ

㻜㻚㻞

0
-1

-0.5

0

-0.5

㻜
㻜㻚㻟

㻜㻚㻠

㻜㻚㻡

㻜㻚㻢

㼆㻌

-1

(b)
Fig. 8Support structure volume, the number of regions
to adhere and the number of subdivision.

Fig. 9 shows the minimal cost value Gmin (N ) withrespect to the number of subdivision N . When theweighting factor α is zero, the cost value monotonously
decreases with the number of subdivision. As theweighting factor α increases, the minimum cost value
is obtained at lower N value.

0.5

1

0
1
2
3
4
5
6
7
8
9
10
11
12 1.5
13
14
15
16
17
18

Sampled deposition direction.

Fig. 10

Deposition direction and cost value.

pled deposit directions. In this experiment, G took theminimum value at direction #16. The optimal numberof subdivision planes was 2. Fig. 11 shows the subdivided objects. The processing time for 19 directions
was 9.2 second on a Windows PC (Intel Core i7-3770K@3.50GHz, GeForce GTX 690). Table 1 shows themeasured processing time for one direction. It took a
few hours to print the objects and we think our optimization process is fast enough.

㻜㻚㻢
㻜㻚㻡

㻯㼛㼟㼠㻌㼢㼍㼘㼡㼑㻌㻳㻌

α =0
㻜㻚㻠

Table 1 PROCESSING TIME
Calculation of
Calculation of
Calculation of
depth map hs(z)Hs(Z,Z’) and M(Z)DP-table

α =0.005

0.18 s㻜㻚㻟

0.27 s<0.001 s

α =0.01㻜㻚㻞

α =0.015

㻜㻚㻝
㻜
㻝

㻞

㻟

㻠

㻡

㻢

㻣

㻤

㻺㻌
Fig. 9The number of subdivision plane and cost value.

Fig. 10-a shows G values with respect to each samPaper

5. 3 Examples
We created physical objects from the Stanford Buddha, Dragon, and Armadillo models12) using the proposed optimization method. The created objects areshown in Fig. 12, Fig. 13 and Fig. 14. In the experiments, we used an FDM 3D printer, named “UP!
Plus 3D Printer”13) . We set the parameter α = 0.01 in
this experiment. As shown in Fig. 12-a, Fig. 13-a andFig. 14-a, when objects were not subdivided, support(5) 5

(a) 1-plane subdivision
Fig. 11(b) 2-plane subdivision (optimal)

Divided 3D CG object.

Fig. 14

Experimental results of the Stanford Armadillo.

Table 2SUPPORT REMOVE TIME AND ADHERING TIME.

BuddhaDragonArmadilloFig. 12Experimental results of the Stanford Buddha.

Fig. 13

Experimental results of the Stanford Dragon.

ing structures could not be completely removed and theobject surface remained rough. In Fig. 12-b, Fig. 13-b
and Fig. 14-b, on the other hand, the optimization resulted in smoother object surfaces and improved qualityof object finish.
6 (6)Number of Support
subdivsion remove(min)

Adhering
(min)

Print
time
(min)

0
19
2(optimal) 90
191(optimal) 1
0
33(optimal) 2

0
5
0
1
0
13

210
215
70
72
127
111

Table 2 shows printing time and working time required in removing supports and adhering subdivided
parts. The table indicates that the optimization largely
decreased support removing time. Subdivision introduces adhering time, which was rather modest with the
Dragon and the Buddha. However, the adhering timewith the Armadillo was relatively large because the subdivided model contained many small regions to adhere.
It was suggested to refine the cost function to avoid
such small subdivided area. The differences in printing
time were not significant.
5. 4 Discussion
Through the experiments, we noticed that the costfunction defined by Eq. (8) may not necessarily be
the best measure of quality and working time but thatthe area of adhering regions and the length of support
structures can be important factors as well. From these
considerations, we have examined some other cost functions based on region areas, support length, and their
combination. Although we have made several pilot experiments using different cost functions and various balancing values (α), we cannot obtained a significant im-

ITE Transactions on Media Technology and Applications Vol. xx, No. xx (20xx)

provement in terms of working time and quality of finish
so far. We realized that more systematic investigations
are necessary to refine the cost function.
The proposed methods digitizes the depth of dividing
planes and the deposit directions in a fixed number oflevels. With respect to the direction, the cost function
cannot be represented by an analytic function in general and it is hard to handle it in a continuous space.
However, it is possible to apply a coarse-to-fine subdivision approach to obtain a good trade-off between the
resolution and cost. We can further subdivide the positions/directions near the minimum values at a higher
sampling rate so that we can search better values in a
depth first fashion.
The proposed method adopts the depth peelingmethod and dynamic programing. To make comparisons with more naive methods, we first implemented
a scan-conversion program that creates sorted lists ofdepth values from which a multi-layered depth image ismade. Surprisingly, the execution time to build a multilayered depth image of the Buddha model at 64 × 64
resolution was 0.031 sec, much smaller than that for thedepth peeling method, 0.18 sec. However, we would still
prefer the depth peeling method to the software scanconversion because of simplicity of implementation. We
also implemented an optimization code based on an exhaustive search. The execution time rapidly grows up
with respect to the maximum number of dividing planes
N and was 13 seconds and 119 seconds for N = 7, 8,
respectively.

6. Conclusion
Although FDM 3D printers have been getting popular, it is still difficult for novice users to determine basic parameters such as deposition directions and model
subdivision. This paper describes an optimization
method that determines the optimal model subdivisionand deposition direction by utilizing CG techniques andthe dynamic programing. The cost function we used isdefined by a sum of supporting structure volumes andthe number of regions to adhere.
Supporting structure volumes can be efficiently computed from multi-layered depth maps on GPU. The optimization can be efficiently achieved by the dynamic
programing. Using the optimization method, we created physical objects with an FDM 3D printer. The
optimization resulted in smoother object surfaces andimproved quality of object finish. It also reduced working time of support removal. The experiments also sugPaper

gested that more refinement of the cost function is necessary to reduce adhering time more effectively. Future
work includes systematic studies on the cost functions.
References
1) ZHU, L., XU, W., SNYDER, J., LIU, Y., WANG, G., AND GUO,
B. 2012. Motion-guided mechanical toy modeling. ACM Trans.
Graph. 31, 6, 127.
2) BICKEL, B., B¨ACHER, M., OTADUY, M. A., LEE, H. R., PFISTER, H., GROSS, M., AND MATUSIK, W. 2010. Design andfabrication of materials with desired deformation behavior. ACM
Trans. Graph. 29, 4.
3) LI, H., VOUGA, E., GUDYM, A., LUO, L., BARRON, J. T.,
AND GUSEV, G. 2013. 3D self-portraits. ACM Trans. Graph. 32,
6 (November), 187.
4) Q. Zhou, J. Panetta, D. Zorin, Worst-case structural analysis,
ACM Trans. on Graphics, vol. 32, No.4, 137, 2013.
5) WANG, W., WANG, T. Y., YANG, Z., LIU, L., TONG, X.,
TONG, W., DENG, J., CHEN, F., AND LIU, X. 2013. Costeffective printing of 3D objects with skin-frame structures. ACM
Trans. Graph. 32, 6, 177:1-177:10.
6) Prévost, R., Whiting, E., Lefebvre, S., and Sorkinehornung, O.
2013. Make it stand: balancing shapes for 3D fabrication. ACM
Trans. Graph. 32, 4 (July), 81:1-81:10.
7) HUANG, X., YE, C., MO, J., AND LIU, H. 2009. Slice data based
support generation algorithm for fused deposition modeling. Tsinghua Science and Technology 14, S1, 223-228.
8) HEIDE, E., 2011. Method for generating and building support
structures with deposition-based digital manufacturing systems,
07. US Patent 20110178621 A1.
9) Z. Hu, K. Lee Concave edge-based decomposition for hybrid rapid
prototyping, International Journal of Machine Tools and Manufacture, vol. 45, No. 1, pp. 35-42,1, January 2005.
10) C. Everitt, Interactive order-independent transparency, Technical report, NVIDIA Corporation, 2001.
11) CGAL, Computational Geometry Algorithms Library,
http://www.cgal.org
12) The Stanford 3D Scanning Repository, http://graphics.stanford.edu/data/3Dscanrep/13) UP! Plus 3D Printer, http://3dprintingsystems.com/products/up-plus-3d-printer/

Satoshi Nakajima received a B.S. and an
M.S degree from Toho University in 2012 and 2014.
He currently works with Fujitsu Broad Solution &Consulting, Inc.

Mikio Shinya is currently a Professor in the
Department of Information Science, Toho University. He received a BSc in 1979, a MS in 1981, anda PhD in 1990 from Waseda University. He joined
NTT Laboratories in 1981, and moved to Toho University in 2001. He was a visiting scientist at the
University of Toronto in 1988-1989. His researchinterests include computer graphics and visual science.

Takashi Yuasa received a B.S. degree fromToho University in 2014. He is a master student at
Graduate School of Science, Toho University.

Michio Shiraishi received the B.A., M.A.,
and Ph.D. degrees from The University of Tokyo in
1997, 1999, 2003, respectively. He is an associate
professor in the Department of Information Science,
Toho University. His current research interests are
non-photorealistic rendering and visualization.

(7) 7