Modelos Massa-Mola para Visualização de Tecidos Flexı́veis
J. M. Pilato Jr., J. S. Espinoza Ortiz, G. A. Giraldi
Laboratório de Visualização Cientı́fica e Realidade Virtual
Coordenação de Ciência da Computação
Laboratório Nacional de Computação Cientı́fica
Caixa Postal 95.113 - 25651-070 – Petrópolis – RJ – Brasil
{pilato,jsespino,gilson}@lncc.br
Abstract. Visual models of flexivel objects have received special attention from
the scientific visualization community. Biologic soft tissue could be considered
a class of flexivel objects, then modeled making use of discret model technics;
too. In this work we present a physical model awhich is also visual implemented
(in real time) where we simulate a flexivel object coliding against a solid one,
by using the spring-mass type discret model.
Resumo. A modelagem visual de objetos flexı́veis tem recebido atenção especial por parte da comunidade de visualização cientı́fica. Por sua vez, os tecidos biológicos podem ser considerados uma classe de objetos flexı́veis, e como
tais também modelados fazendo uso de técnicas de modelagem discreta. Neste
artigo apresentamos um modelo fı́sico e a sua implementação visual (em tempo
real) onde simulamos a colisão de um objeto flexı́vel e um objeto sólido, fazendo
uso do modelo discreto do tipo massa-mola.
1. Introdução
O realismo é um elemento fundamental quando trata-se de visualização de fenômenos
fı́sicos em computação gráfica. O nosso interesse está focado na modelagem dos tecidos
flexı́veis do tipo biológico. É com essa finalidade que procuramos por modelos realistas,
mas que também representem as propriedades fı́sicas do objeto em questão.
No caso dos tecidos biológicos, as suas propriedades mecânicas dependem da sua
composição mineral. É notório que em função do seu conteúdo mineral podemos em
geral distinguir duas classes de tecidos biológicos. Por exemplo, ossos e dentes contêm
minerais e formam parte do grupo que é denominado de tecidos duros. Enquanto que,
pele, músculo, artérias e pulmões formam um segundo grupo, chamado de tecidos suaves.
Estes últimos não contêm minerais e por esta razão são bem mais deformáveis do que os
tecidos duros.
Tem sido bastante comum modelar o comportamento mecânico dos tecidos suaves
aplicando técnicas usadas na modelagem de materiais eláticos do tipo borracha. Neste
caso, teorias de deformação finita são freqüentemente usadas para descrever seu comportamento mecânico. Entretanto, existem diferenças marcantes no que diz respeito à estrutura material entre os tecidos suaves e os materiais elásticos, principalmente no modo
em que eles respondem sob tensão aplicada. Os tecidos suaves facilmente alcançam
consideráveis extensões sob um nı́vel de tensão relativamente baixo para após tornarse mais rı́gidos sob tensões de nı́veis mais intensos. Esse comportamento é mostrado
na Figura (1), onde são comparados as diferentes respostas à tensão aplicada a materiais
elásticos como a borracha (figura à esquerda) com a do tecido suave (figura à direita). Por
outro lado, a borracha é tipicamente um material isotrópico, em contraste com os tecidos
suaves que manifestam uma tı́pica anisotropia por causa da sua composição de fibras de
Colágeno[Ogden 2001] [Holzapfel 2000][Holzapfel 2001].
Figura 1. Tı́pica resposta da deformação a uma simples tensão sob a borracha
(a) e sob o tecido suave (b). Note que a tensão nominal (t) é positiva, por tanto
a deformação λ é também positiva, sendo esta última definida em unidades do
comprimento original da amostra.
Aplicações em simulações cirúrgicas tem aumentado bastante a demanda em
pesquisa de modelagem do comportamento mecânico de tecidos biológicos suaves, procurando por um cálculo preciso e em tempo real da deformação do tecido. Uma das
técnicas amplamente utilizadas está baseada em modelos constitutivos que toma conta
da natureza contı́nua dos tecidos, onde técnicas como as de elementos finitos freqüentemente auxiliam para uma correta simulação do comportamento mecânico dos tecidos
[Bro-Nielsen 1998][Keeve et al. 1996][Koch et al. 1996]. Por outro lado, como alternativa, também tem sido utilizados os modelos discretos baseados em modelos massa-mola
[Mosegaard 2003][Teschrer et al. 2000]. Pesquisas mostram que métodos de Elementos
Finitos baseados na aproximação do contı́nuo são apropriados para descrever de modo
mais realı́stico o comportamento mecânico dos tecidos suaves, entretanto, modelos baseados em massa-mola são computacionalmente econômicos, se comparados com os anteriores.
Neste artigo apresentamos um modelo fı́sico discreto que utiliza sistemas de
massa-mola. A Seção 2 é iniciada como uma apresentação simples do esquema massamola comunmente usado em visualização cientı́fica, com a finalidade de representar tecidos (objetos) flexı́veis inorgânicos. Nestes modelos são geralmente considerados efeitos
lineares, tanto como para a elasticidade da mola assim como para o atrito entre massas. Imediatamente é apresentada de forma breve uma extensão da teoria onde podem
ser incluı́das interações do tipo viscoelástico assim como interações elásticas do tipo nãolinear, que podem ser úteis na modelagem de tecidos biológicos. Já na Seção 3 nos especializamos nos modelos lineares onde consideramos um sistema composto por massas
conectadas através de molas representando um objeto flexı́vel, esquematizado de forma
que irá colidir com um obstáculo sólido fixo. Diferentes situações são consideradas. Logo
na Seção 4 são apresentados detalhes da implementação visual dos nossos experimentos.
Finalmente na Seção 5 apresentamos as nossas conclusões assim como as perspectivas
futuras.
2. O modelo massa-mola
Consideremos um corpo fı́sico flexı́vel que é composto por um conjunto de partı́culas de
massas mi conectadas através de molas. Estabelece-se também que inicialmente ambas,
as forças elásticas e de atrito entre partı́culas, reagem linearmente. Por esta razão, vamos considerar que o corpo elástico obedece e lei de Hooke. Consideremos agora que o
mesmo corpo seja suportado a um sistema de referência inercial determinado e sujeitado
a um conjunto de forças Qi que age nos pontos i ( = 1, 2, 3, . . . , n). Em seguida consideramos uma deflexão generalizada qi no ponto i, em direção a força Qi , sendo linearmente
proporcional as forças Q1 , Q2 , Q3 , . . . , Qn , e vice versa, para i ∈ (1, 2, 3, ..., n) teremos:
qi = Ci1 Q1 + Ci2 Q2 + . . . + Cin Qn ,
Qi = Ki1 q1 + Ki2 q2 + . . . + Kin qn .
(1)
Onde as constantes de proporcionalidade Cij , Kij são independentes das forças
e deslocamentos. Cij e Kij são os coeficientes de flexibilidade e de rigidez, respectivamente. O significado fı́sico de Kij é a força requerida para agir no ponto i devido a um
deslocamento unitário no ponto j, quando uns outros pontos (diferentes do que j) são
fixos. Analogamente, Cij é a deflexão em i devido a uma força unitária que age em j. As
equações acopladas, acima, implicam que existe um estado natural de equilı́brio ao qual
o corpo retornará sempre que todo as forças externas forem removidas, e onde o princı́pio
do superposição de forças se aplica. Além disso, implica que o trabalho total feito por um
conjunto de forças não depende da ordem em que as forças são aplicadas. Este trabalho
feito por cada força é 21 Qi qi , de maneira que o trabalho total feito pelo conjunto, que é
armazenado como energia de deformação no corpo elástico, é dado por:
W =
n
n X
n
n X
n
1X
1X
1X
Qi qi =
Kij qi qj =
Cij Qi Qj .
2 i=1
2 i=1 j=1
2 i=1 j=1
(2)
As matrizes Kij e Cij são simétricas. Em outras palavras, o deslocamento em um
ponto i devido a uma força unitária que age em um outro ponto j são iguais ao deslocamento em j devido a uma força unitária que age em i, isto é, são positivos e na mesma
direção em cada ponto.
Agora podemos proceder a escrever a equação do movimento de um grupo de
partı́culas de massas mi contidas em um corpo elástico. O princı́pio de D’Alambert estabelece que as partı́culas podem ser consideradas em seu estado de equilı́brio se as forças
inerciais forem aplicadas no sentido invertido sob as partı́culas. Assim, se o corpo estiver
em estado estacionário e nós considerarmos forças de atrito lineares, então a equação
dinâmica escrita na forma vetorial é
mi q̈i +
n
X
j=1
Kij qj +
n
X
βij q̇j − Fie = 0; (i = 1, 2, ..., n).
(3)
j=1
Aqui o q̈i é a aceleração resultante da i−ésima partı́cula, em seguida estão sendo
considerados as forças agindo sob a partı́cula em questão: como as forças internas do tipo
elástico (equação 1), assim como outras forças internas de amortecimento, onde qj , q̇j são
respectivamente o deslocamento e a velocidade da j−ésima partı́cula, e Fie é uma força
inercial externa.
Entretanto, a força de amortecimento poderia ser muito mais complexa do que
a força de atrito acima equacionada, por exemplo pode ser aerodinâmica na origem, ou
com viscoelasticidade não-linear, tal e como acontece com os tecidos suaves. Para uma
generalização imediata da teoria, vamos fazer uso da mais geral das equações constitutivas utilizada na modelagem de sólidos viscoelásticos, onde podemos descrever as forças
internas agindo sobre uma determinada partı́cula, localizada num certo ponto i, por meio
de um deslocamento qj localizado num ponto j, como
Fij (t) =
Z
t
−∞
Gij (t − τ )
dqj (τ )
dτ .
dτ
(4)
Seguindo a formulação de Boltzman [Fung 1994], Gij (t) é chamada de função de
relaxação, expressada matematicamente da seguinte forma:
Gij (t) = Kij + βij e−tν1 + ... + γij e−tνn ,
(5)
onde νi são chamadas de freqüências de relaxação, Kij , βij , γij são constantes.
Desta forma, considerando tanto as forças do tipo interna quanto as externas que agem
sob a partı́cula, a Equação (3) é re-escrita da seguinte maneira:
mi
d2 qi (t)
dt2
+
Rt
j=1 −∞
Pn
Gij (t − τ )
dqj (τ )
dτ
dτ
(e)
− Fi
= 0.
(6)
3. Implementação Numérica do Modelo
A presente implementação gráfica consiste de um sistema massa-mola, inicialmente a
uma altura pré-determinada h a respeito de uma superfı́cie referencial onde encontrase um objeto sólido em repouso. Instantes depois o objeto flexı́vel irá colidir com o
objeto sólido por causa da força da gravidade, conforme ilustra a Figura (2). Além das
interações elásticas lineares massa-mola são também consideradas forças de atrito do tipo
Newtoniana, como as da equação (3), porém proporcionais à velocidade instantânea da
respectiva partı́cula.
Em seguida nos focamos na i−ésima partı́cula alocada no ponto i da malha. A
equação dinâmica para a presente partı́cula será determinada pela Equação (3) onde: o
primeiro termo a esquerda é a sua aceleração resultante, o segundo termo representa
a força elástica do tipo Hooke (onde Kij é a constante de elasticidade), o seguinte
termo constitue a força de atrito, Fie é a força gravitacional (mi g). Imediatamente as
equações diferenciais de segunda ordem são resolvidas fazendo uso do esquema denominado Leap-Frog, que é um método iterativo de segunda ordem[Müller et al. 2003], o qual
consiste em: Dados os deslocamentos e velocidades para um determinado instante “ t”, a
aceleração respectiva é determinada a partir da equação (3), em seguida os novos deslocamentos e velocidades para um novo instante de tempo “ t + ∆ t” serão computados por
meio das seguintes relações,
Figura 2. Esquema objeto flexı́vel - objeto sólido, rearranjado para uma colisão.
h e g são a altura e a aceleração da gravidade, respectivamente.
~q˙ i (t + ∆ t) = ~q˙ (t)i + ∆2 t ¨~qi (t) ;
~qi (t + ∆ t) = ~qi (t) + ∆ t~q˙ i (t + ∆ t) .
(7)
A seguir apresentamos os resultados de nossos experimentos em três situações
distintas. Na primeira situação o objeto flexı́vel é composto por um conjunto de massas homogêneas; 1681 em número, onde estabelecemos os seguintes parâmetros: mi =
0.5, K = 130, βij = 0.8, ∆ t = 10−2 , ~g = −10~j e q̇i (0) = −10~j, sendo ~j = (0, 1, 0), de
maneira que o objeto flexı́vel foi lançado contra o objeto sólido, representado por meio
de uma esfera fixa no espaço.
Na Figura (4) de (a)-(d) mostra-se a seqüência da queda do objeto flexı́vel, onde
apenas é mostrada a deformação da malha. Já na seqüencia (e)-(h), da mesma figura,
visualizamos a mesma situação anterior mas destacando a deformação das molas, segundo
a escala de cores abaixo. Note que da esquerda à direita as deformações incrementam-se
seguindo as mudanças de cores.
Figura 3. Escala de deformação: As cores extremas da esquerda e direita representam relaxação e deformação máxima, respectivamente.
A seguinte situação apresentada, Figura (5), difere da anterior no que diz respeito
à homogeneidade das massas, pois neste caso usamos valores diferentes para as mesmas.
Para isso dividimos o objeto flexı́vel em duas regiões simétricas, onde em uma delas a
massa é o dobro da outra, ou seja, mamarelo
= 0.5 e mverde
= 1.0, respectivamente.
i
i
Exceto as massas, os parâmetros restantes serão os mesmos da primeira situação. Neste
caso o experimento também mostra ambas as deformações da malha (à esquerda) assim
como as deformações seguindo a escala de cores acima mencionada (à direita).
Um último experimento apresenta o objeto flexı́vel, novamente considerando as
massas homogêneas, mas desta vez colidindo com dois objetos sólidos posicionados nas
duas extremidades do objeto flexı́vel. Uma vez que o objeto flexı́vel se choca com os objetos sólidos este passa a se deslizar entre as esferas para naturalmente acabar se soltando
do obstáculo. A presente situação pode ser observada na Figura (6), desta vez apenas
(a)
(e)
(b)
(f)
(c)
(g)
(d)
(h)
Figura 4. Na sequência instantes da queda do objeto flexı́vel em (a)-(d). A mesma
sequência é mostrada em (e)-(h) onde as deformações seguem a escala de cores
da Figura (3).
(a)
(e)
(b)
(f)
(c)
(g)
(e)
(h)
Figura 5. Um outro experimento similar ao mostrado na figura anterior, onde as
cores diferenciadas representam massas com valores distintos. A seqüência (a)(e) mostra as deformações da malha, enquanto que a seqüência (d)-(f) mostra as
mesmas seguindo a escala de cores.
visualizamos as deformações segundo a escala de cores1 .
4. Implementação Visual
O objeto flexı́vel consiste em uma malha estruturada de 1681 pontos, numa configuração
matricial de 41x41 conectadas conforme ilustra a Figura (7). Por sua vez, o corpo rı́gido
pelo qual o objeto elástico se choca é representado por uma esfera, ou um par delas.
Para cada instante de tempo da simulação calculamos: a aceleração, a velocidade e
a nova posição de cada massa, conforme as equações (3) e (7). Com base nestas equações
utilizamos o algoritmo (1) para simular a cada instante de tempo a dinâmica do objeto
flexı́vel.
Algoritmo 1 : Update()
para cada instante faça
para cada Pi,j faça
Calcular aceleração de Pi,j ;
fim para
para cada Pi,j faça
Calcular velocidade de Pi,j ;
Calcular nova posição de Pi,j ;
fim para
Renderizar o modelo;
fim para
Para desenvolver a simulação fizemos uso da ferramenta Visualization Toolkit
(VTK) [Schroeder et al. 1997]. O VTK é uma biblioteca livre, de código aberto, amplamente utilizado em centros de pesquisa com a finalidade de desenvolver aplicações em
computação gráfica, processamento de imagens e visualização cientı́fica. O VTK consiste
em uma biblioteca de classes em C++, cujo projeto de desenvolvimento e execução são
fortemente influenciados por princı́pios de orientação a objetos.
A fim de manter melhor compatibilidade com o VTK, e aproveitar das inúmeras
vantagens do paradigma de orientação a objetos [Silva Neto et al. 2005], desenvolvemos
nossa simulação em C++. O Diagrama de classes da Figura (8) representa, de forma
estruturada, a implementação da simulação realizada nesse trabalho.
No diagrama de Classes, MassSpring é a principal classe da simulação, foi nela
que implementamos o algoritimo (1). Para isso, a mesma possui os membros MassCollection e SpringCollection, que são uma lista de massas (classe Mass) e molas (classe
Spring), respectivamente. Naturalmente, a classe Mass representa cada massa presente
no sistema, possuindo atributos como: coordenada, velocidade e aceleração. Por sua vez,
a classe Spring representa as molas e é usada para identificar a conectividade entre as
partı́culas.
5. Conclusões
O esquema de colisão objeto flexı́vel-objeto sólido mostrado no presente artigo apresenta
um nı́vel realı́stico razoável. As deformações da malha na escala de cores mostram clara1
Imagens dos experimentos, aqui mostrados, em tempo real poderão ser disponibilizadas aos leitores.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figura 6. Deslizamento do objeto flexı́vel entre as esferas, apenas são mostradas
as deformações da malha segundo a escala de cores. Neste caso os parâmetros
são os mesmos considerados na Figura (4).
Figura 7. Esquema de conexão da malha utilizada na simulação.
Figura 8. Diagrama de Classes da simulação de tecido.
mente maiores deformações na parte superior da esfera, como no caso das Figuras (4-h) e
(5-h), sugerindo que aquela região sofre elevadas tensões por conta do sistema restante. Já
no caso da Figura (6), esta indica uma seqüência de deslizamento fisicamente apropriada.
No entanto consideramos que há ainda a possibilidade de fazermos melhoras na
nossa metodologia, no que diz respeito ao cálculo numérico. O esquema iterativo utilizado
é razoável apenas para tempos de cálculo não muito longos. Nesse caso poderı́amos usar
esquemas de cálculo mais elaborados, como do tipo Runge Kutta [Press et al. 1992], testando então diferentes ordens de aproximação à procura de precisão assim como também
por eficiência. Um outro ponto crucial a respeito da modelagem e visualização de tecidos biológicos consiste na utilização da formulação que faz uso da Equação (6), onde
deveremos resolver equações acopladas integro-diferenciais do tipo Volterra, as quais já
estamos implementando. É desta forma que pretendemos continuar nossa linha de trabalho visando por uma visualização ótima e realı́stica do comportamento mecânico dos
tecidos biológicos.
Agradecimentos
Agradecemos ao Conselho Nacional de Desenvolvimento Cientı́fico e Tecnológico CNPq, pelo apoio financeiro dado a este projeto.
Referências
Bro-Nielsen, M. (1998). Finite element modelling in surgery simulation. Proc. of the
IEEE: Special Issue on Virtual and Augmented Reality in Medicine 86 (3) 524-30.
Fung, Y. (1994). A First Course in Continuum Mechanics. Prentice Hall, Englewood
Cliff, New Jersey, USA.
Holzapfel, G. (2000). Nonlinear Solid Mechanics. Chichester: Wiley.
Holzapfel, G. A. (2001). Biomechanics of Soft Tissue. In Lecture, J.ed., Handbook of
Material Behavior Models. Boston: Academic Press 1049-63.
Keeve, E., Girod, S., Pfeifel, P., and Girod, B. (1996). Anatomy-based facial tissue modelling using the finite element method. Proc. IEEE Visualization, San Francisco, CA.
USA.
Koch, R., Gross, M., Bueren, D., Franhauser, G., Parish, Y., and Carls, F. (1996). Simulation facial surgery using finite element modles. SIGGRAPH’96, ACM Computer
Graphycs, August 30.
Müller, M., Charypar, D., and Gross., M. (2003). Particle-based fluid simulation for
interactive applications. In Proceedings of ACM SIGGRAPH symposium on Computer
animation.
Mosegaard, J. (2003). Lr-spring-mass model for cardiac surgical simulation. Medicine
Meets Virtual Reality, (12) 256-58.
Ogden, R. (2001). Nonlinear Elasticity, Anisotropy, Material Stability and Residual Stress
in Soft Tissue.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1992). Numerical
Recipes in C: The Art of Scientific Computing. Cambridge University Press, New York,
NY, USA.
Schroeder, W., Martin, K., and Lorensen, B. (1997). The Visualization Toolkit: An ObjectOriented Approach to 3D Graphics. Prentice Hall.
Silva Neto, A. A., Rodrigues, P., Giraldi, G. A., and Apolinario Jr., A. L. (2005).
Animação computacional de fluidos via smoothed particle hydrodynamics. Technical Report 06/2005, Laboratório Nacional de Computação Cientı́fica - Petrópolis/RJ.
Teschrer, M., Girod, S., and Girod, B. (2000). Direct computation of nolinear soft-tissue
deformation. Vison Modeling and Visualization VMV’000, Saarbucken, Germany,
November 22-24.