Update
Updated by April 14.
Intermediate result:
CPU part.
1 Surface Mass-Spring Modeling
2 Internal Skeleton Modeling
3 The Hybrid Model
GPU part
4 basic shadering effect.
5 CUDA installation and configuration
6 kernel of the computing
1 Basic Mass-Spring Modeling.
The basic idea of Mass-Spring is modeling a soft body with particles linked by springs. at each iteration, we update the position, velocity and acceleration for each particle based on this equation:
Intermediate result:
CPU part.
1 Surface Mass-Spring Modeling
2 Internal Skeleton Modeling
3 The Hybrid Model
GPU part
4 basic shadering effect.
5 CUDA installation and configuration
6 kernel of the computing
1 Basic Mass-Spring Modeling.
The basic idea of Mass-Spring is modeling a soft body with particles linked by springs. at each iteration, we update the position, velocity and acceleration for each particle based on this equation:
we can write the C++ code like this:
for(i=0;i<(int)springs.size();i++)
{
glm::vec3 p1 = vec3(X[springs[i].p1]);
glm::vec3 p1Last = vec3(X_last[springs[i].p1]);
glm::vec3 p2 = vec3(X[springs[i].p2]);
glm::vec3 p2Last = vec3(X_last[springs[i].p2]);
glm::vec3 deltaP = p1-p2;
float dist = glm::length(deltaP);
float leftTerm = -springs[i].Ks * (dist-springs[i].rest_length);
glm::vec3 springForce = leftTerm *glm::normalize(deltaP);
F[springs[i].p1] += springForce;
F[springs[i].p2] -= springForce;
}
In the above code for each springs, we update the two particles linked by this spring. we use the Euler method to update the positions for each particle.
for(i=0;i<(int)springs.size();i++)
{
glm::vec3 p1 = vec3(X[springs[i].p1]);
glm::vec3 p1Last = vec3(X_last[springs[i].p1]);
glm::vec3 p2 = vec3(X[springs[i].p2]);
glm::vec3 p2Last = vec3(X_last[springs[i].p2]);
glm::vec3 deltaP = p1-p2;
float dist = glm::length(deltaP);
float leftTerm = -springs[i].Ks * (dist-springs[i].rest_length);
glm::vec3 springForce = leftTerm *glm::normalize(deltaP);
F[springs[i].p1] += springForce;
F[springs[i].p2] -= springForce;
}
In the above code for each springs, we update the two particles linked by this spring. we use the Euler method to update the positions for each particle.
2 Soft tissue modeling and Internal Skeleton Modeling
Since soft tissue is a special mass-spring model, during the implementation, we improve the structure of traditional Mass-Spring system in order to refine the effect of simulation. Firstly, we add a curvature force. Curvature force controls the degree of bending and twisting of soft tissue. We imitate the force by bringing in an assistant spring: angular spring. As shown below P1, P2 and P3 are the surface points of soft tissue. The angular spring links point A with the mid point of P2P3. Secondly, we induct a concept of fixed position in order to prevent soft tissue from escaping from the original location. Here, we presume that each point has a corresponding fixed spot located in the original position of the point, and the point connects with the spot through springs named return-springs. As a result, soft tissue always has a tendency to go back to the original position.
Since soft tissue is a special mass-spring model, during the implementation, we improve the structure of traditional Mass-Spring system in order to refine the effect of simulation. Firstly, we add a curvature force. Curvature force controls the degree of bending and twisting of soft tissue. We imitate the force by bringing in an assistant spring: angular spring. As shown below P1, P2 and P3 are the surface points of soft tissue. The angular spring links point A with the mid point of P2P3. Secondly, we induct a concept of fixed position in order to prevent soft tissue from escaping from the original location. Here, we presume that each point has a corresponding fixed spot located in the original position of the point, and the point connects with the spot through springs named return-springs. As a result, soft tissue always has a tendency to go back to the original position.
In order to establish the internal model, Distance Mapping Method is employed to calculate the skeleton, along which media atoms are selected evenly. Then we simplify M-Rep with the purpose of reducing the model’s complexity. We reset the topology of medial atom and implied boundary just like hub and spokes.
3 The Hybrid Model
The hybrid model is established after the above processes. Then the skeleton and spokes are considered to be springs. Finally, we should set the parameters of the model appropriately in order to simulate soft tissue effectively and efficiently. we use a simple techniques of finding the closest surface particle to link with the spoke.
3 The Hybrid Model
The hybrid model is established after the above processes. Then the skeleton and spokes are considered to be springs. Finally, we should set the parameters of the model appropriately in order to simulate soft tissue effectively and efficiently. we use a simple techniques of finding the closest surface particle to link with the spoke.
GPU part.
4 basic shadering effect.
We use the fixed pipeline in OpenGL to render the scene. The effect is as above.
5 CUDA installation and configuration
The platform is VS2012 and the CUDA version is 5.0 which is the newest version. The installation and configuration of these two platforms is pretty tricky, since the CUDA 5.0 is not supportive to VS2012. Some codes in .lib and .h have been modified to integrated to VS2012.
6 kernel of the computing
Currently, I'm still working on this part, the basic idea is each thread is responsible for one particle, the neighboring particles are used to update the position of this particle. The up to date code is as follows:
__global__ void verlet( float4 * pos_vbo, float4 * g_pos_in, float4 * g_pos_old_in, float4 * g_pos_out, float4 * g_pos_old_out,
int2 texsize, float2 step, float damp, float mass, float dt, float2 inv_cloth_size,vector<ajacent> ajacent)
{
uint index = blockIdx.x * blockDim.x + threadIdx.x;
volatile float4 posData = g_pos_in[index];
volatile float4 posOldData = g_pos_old_in[index];
float3 pos = make_float3(posData.x, posData.y, posData.z);
float3 pos_old = make_float3(posOldData.x, posOldData.y, posOldData.z);
float3 vel = (pos - pos_old) / dt;
const float3 gravity=make_float3(0.0f,-0.00981f,0.0f);
float3 force = gravity*mass + vel*damp;
for(int i=)
{
//TODO here, using the code above to update the force based on the adjacent particles
}
float3 acc = make_float3(0, 0, 0);
if(mass!=0)
acc = force / mass;
// verlet
float3 tmp = pos;
pos = pos * 2 - pos_old + acc * dt * dt;
pos_old = tmp;
syncthreads();
pos_vbo[index] = make_float4(pos.x, pos.y, pos.z, posData.w);
g_pos_out[index] = make_float4(pos.x, pos.y, pos.z, posData.w);
}
4 basic shadering effect.
We use the fixed pipeline in OpenGL to render the scene. The effect is as above.
5 CUDA installation and configuration
The platform is VS2012 and the CUDA version is 5.0 which is the newest version. The installation and configuration of these two platforms is pretty tricky, since the CUDA 5.0 is not supportive to VS2012. Some codes in .lib and .h have been modified to integrated to VS2012.
6 kernel of the computing
Currently, I'm still working on this part, the basic idea is each thread is responsible for one particle, the neighboring particles are used to update the position of this particle. The up to date code is as follows:
__global__ void verlet( float4 * pos_vbo, float4 * g_pos_in, float4 * g_pos_old_in, float4 * g_pos_out, float4 * g_pos_old_out,
int2 texsize, float2 step, float damp, float mass, float dt, float2 inv_cloth_size,vector<ajacent> ajacent)
{
uint index = blockIdx.x * blockDim.x + threadIdx.x;
volatile float4 posData = g_pos_in[index];
volatile float4 posOldData = g_pos_old_in[index];
float3 pos = make_float3(posData.x, posData.y, posData.z);
float3 pos_old = make_float3(posOldData.x, posOldData.y, posOldData.z);
float3 vel = (pos - pos_old) / dt;
const float3 gravity=make_float3(0.0f,-0.00981f,0.0f);
float3 force = gravity*mass + vel*damp;
for(int i=)
{
//TODO here, using the code above to update the force based on the adjacent particles
}
float3 acc = make_float3(0, 0, 0);
if(mass!=0)
acc = force / mass;
// verlet
float3 tmp = pos;
pos = pos * 2 - pos_old + acc * dt * dt;
pos_old = tmp;
syncthreads();
pos_vbo[index] = make_float4(pos.x, pos.y, pos.z, posData.w);
g_pos_out[index] = make_float4(pos.x, pos.y, pos.z, posData.w);
}
Final Project Proposal
What is the goal of the project?
The project is an implementation of CUDA based soft tissue deformation simulation. The deformation model is spring mass model. To speed up the simulation performance for surgical tissue deformation, CUDA based GPU computing is adopted, while related data structures and algorithm are designed and implemented for the parallel computation.
What is your motivation behind this project?
The medical training systems based on virtual simulation are highly desired since minimally invasive surgical techniques have become popular to patients. The training system helps surgeon trainees to acquire, practice and evaluate their surgical skills, and the key component of such a system is to simulate the dynamic procedure such as 3D biological tissue deformation in surgical operation.
What is the prior work or related work that your work builds on? Basically we want to make sure that you are aware of state of the art.
S. Zhang, L. Gu1, P. Huang, and J. Xu, “Real-time simulation of deformable soft tissue based on mass-spring and medial representation,” in Computer Vision for Biomedical Image Applications: First International Workshop, CVBIA 2005.
X. Provot, “Deformation constraints in a spring-mass model to describe rigid cloth behavior,” in Proceedings of Graphics Interface, 1995, pp. 147–154.
Mass–spring–damper modelling of the human body to study running and hopping – an overview.
4. What do you plan to accomplish during the course? There would be a final project presentation (sometime between April 26 and May 3, i.e. before the final) and you are expected to show the results or give a writeup. So how much can you accomplish in this time frame. In the past, it is not uncommon for the students to underestimate the effort or project goals.
Firstly, I will implement a basic spring mass model, create the rendering effect in order to make it more realistic, after them has been done, I will leverage CUDA to accelerate the simulation. After that, I will analysis the efficiency of CUDA simulation. (No complex collision detection is needed, since we only have to simulate the knife interacted with soft tissue)