척추의 의료 영상 기반 유한요소해석

다음 제품이 필요합니다.

이 예제에서는 FE(유한 요소) 해석을 사용하여 축 방향 압축 하에 있는 척추골의 뼈 응력과 변형률을 추정하는 방법을 보여줍니다.

뼈는 중력, 움직임, 근육 수축으로 인해 기계적 하중을 받습니다. 뼈 조직 내의 응력과 변형률을 추정하면 뼈 강도와 골절 위험을 예측하는 데 도움이 될 수 있습니다. 이 예제에서는 영상 기반 FE 해석을 사용하여 축 방향 압축 하에 있는 척추골 내의 뼈 응력과 변형률을 예측합니다. 영상 기반 FE에서는 CT(컴퓨터 단층촬영) 스캔과 같은 의료 영상을 사용하여 FE 모델의 형상 및 물질 속성을 생성합니다.

이 예제에서는 다음 워크플로를 사용하여 축 방향 하중을 받는 단일 척추골의 3차원 모델을 만들고 분석합니다.

데이터 다운로드하기

이 예제에서는 DICOM 파일 디렉터리로 저장된 CT 스캔을 사용합니다. 이 스캔은 3개의 CT 볼륨을 포함하는 데이터 세트에 속합니다. 다음 코드를 실행하여 MathWorks 웹사이트에서 데이터 세트를 zip 파일로 다운로드한 다음 압축을 풉니다. 데이터 세트의 전체 크기는 약 81MB입니다.

zipFile = matlab.internal.examples.downloadSupportFile("medical","MedicalVolumeDICOMData.zip");
filepath = fileparts(zipFile);
unzip(zipFile,filepath)

데이터 세트를 다운로드하고 압축을 풀면 이 예제에 사용되는 스캔을 dataFolder 디렉터리에서 확인할 수 있습니다.

dataFolder = fullfile(filepath,"MedicalVolumeDICOMData","LungCT01");

CT 스캔에서 척추골 분할하기

이 예제에서는 척추의 한 척추골에 대한 이진 분할 마스크를 사용합니다. 마스크는 의료 영상 레이블 지정기 앱을 사용하여 흉부 CT 스캔에서 척추를 분할하여 생성되었습니다. 앱에서 의료 영상 볼륨을 분할하는 방법에 대한 예는 Label 3-D Medical Image Using Medical Image Labeler 항목을 참조하십시오.

Chest CT scan loaderd in Medical Image Labeler app.

분할 마스크는 이 예제에 지원 파일로 첨부되어 있습니다. 마스크를 medicalVolume 객체로 불러옵니다.

labelVol = medicalVolume("LungCT01.nii.gz");

의료 볼륨 객체의 VolumeGeometry 속성에는 스캔의 환자 좌표계를 정의하는 medicalref3d 객체가 포함되어 있습니다. 마스크의 medicalref3d 객체를 새 변수 R로 추출합니다.

R = labelVol.VolumeGeometry;

척추골 형상을 추출하고 FE 메시 생성하기

척추골 마스크의 등가곡면 추출하기

extractIsosurface 함수를 사용하여 척추골 마스크의 등가곡면 면과 꼭짓점을 계산합니다. vertices 좌표는 행, 열, 슬라이스로 정의되는 내재적 좌표계에 있습니다.

isovalue = 1;
[faces,vertices] = extractIsosurface(labelVol.Voxels,isovalue);

intrinsicToWorld 객체 함수를 사용하여 vertices를 R로 정의된 환자 좌표계로 변환합니다. X, Y 및 Z 출력값에 따라 곡면 점의 xyz 좌표가 밀리미터 단위로 정의됩니다.

I = vertices(:,1);
J = vertices(:,2);
K = vertices(:,3);

[X,Y,Z] = intrinsicToWorld(R,I,J,K);
verticesPatient = [X Y Z];

꼭짓점 좌표를 미터로 변환합니다.

verticesPatientMeters = verticesPatient.*10^-3;

척추골 등가곡면을 삼각분할된 Surface 객체로 표시합니다.

triangul = triangulation(double(faces),double(verticesPatientMeters));
viewer = viewer3d;
surface = images.ui.graphics.Surface(viewer,data=triangul,Color=[1 0 0],Alpha=1);

PDE 모델 형상 만들기

3개의 방정식으로 구성된 시스템에 대한 일반적인 PDEModel 모델 컨테이너를 만듭니다. 일반적인 모델을 사용하면 골밀도에 따라 공간적으로 변화하는 물질 속성을 할당할 수 있습니다.

model = createpde(3);

환자 좌표의 등가곡면을 입력값으로 받는 geometryFromMesh (Partial Differential Equation Toolbox) 객체 함수를 사용하여 모델의 형상을 지정합니다.

geometryFromMesh(model,verticesPatientMeters',faces');
pdegplot(model,FaceLabels="on")

Figure contains an axes object. The axes object contains 6 objects of type quiver, text, patch, line.

모델 형상에서 메시 생성하기

FE 메시 요소의 최대 가장자리 길이를 지정합니다. 이 예제에서는 최대 가장자리 길이로 1.6mm를 사용하며, 이는 메쉬 수렴 해석을 통해 결정된 값입니다.

hMax = 0.0016;

모델 형상에서 메시를 생성합니다. 최대 가장자리 길이를 지정하고, 최소 가장자리 길이가 최대 가장자리 길이의 절반이 되도록 지정합니다. 요소는 2차 사면체(quadratic tetrahedra)입니다.

msh = generateMesh(model,Hmax=hMax,Hmin=hMax/2);
pdemesh(model);

Figure contains an axes object. The hidden axes object contains 5 objects of type quiver, text, patch.

물질 속성 할당하기

뼈 물질 속성은 공간적 위치와 방향에 따라 달라집니다.

공간 위치 — 골밀도가 더 높은 국소 영역은 골밀도가 낮은 영역보다 강성이 더 큽니다.
방향 — 뼈의 강성은 하중 방향에 따라 달라집니다. 이 예제에서는 횡방향 등방성을 사용하여 뼈를 모델링합니다. 횡방향 등방성 물질은 횡단면에 비해 축 방향으로 강성이 더 크며 횡단면 내에서는 균일합니다.

이 예제에서는 공간적으로 변화하는 횡방향 등방성 선형 탄성 물질 속성을 할당합니다.

척추골 마스크 내 HU 값 추출하기

공간적으로 변화하는 물질 속성을 할당하려면 CT 강도 값을 FE 메시 좌표에 매핑해야 합니다.

medicalVolume 객체를 사용하여 분할되지 않은 원본 CT 스캔을 불러옵니다. medicalVolume 객체는 자동으로 강도를 하운스필드 단위(HU)로 변환합니다. HU 값은 방사선 밀도를 측정하는 기준입니다.

DCMData = medicalVolume(dataFolder);

척추골 마스크 내의 CT 복셀의 인덱스와 강도를 추출합니다.

trueIdx = find(labelVol.Voxels==1);
HUVertebra = double(DCMData.Voxels(trueIdx));

trueIdx의 선형 인덱스를 복셀 단위의 첨자 인덱스로 변환합니다.

[row,col,slice] = ind2sub(size(labelVol.Voxels),trueIdx);

첨자 인덱스를 환자 좌표로 변환하고 좌표를 밀리미터에서 미터로 변환합니다.

[X2,Y2,Z2] = intrinsicToWorld(R,col,row,slice);
HUVertebraMeters = [X2 Y2 Z2].*10^-3;

CT 복셀에서 메시 노드로 HU 값 매핑하기

FE 노드 좌표는 분산되어 있으며 CT 복셀 그리드와 정렬되어 있지 않습니다. 복셀 그리드와 FE 노드 사이의 3차원 보간을 정의하는 scatteredInterpolant 객체를 만듭니다.

F = scatteredInterpolant(HUVertebraMeters,HUVertebra);

계수 지정하기

모델에 대한 PDE 계수를 지정합니다. 이 구조 해석에서 c 계수는 모델의 물질 강성을 정의하며, 이는 컴플라이언스와 반비례합니다. Partial Differential Equation Toolbox™에서 공간적으로 변화하는 c 계수를 정의하려면 c를 함수 형태로 표현합니다. 이 예제에서 HU2TransverseIsotropy 헬퍼 함수는 골밀도에 따라 횡방향 등방성 물질 속성을 정의합니다. 골밀도는 산점 보간 함수 F를 사용하여 주어진 위치에 대해 계산됩니다. 헬퍼 함수는 익명 함수 ccoeffunc에 래핑되어 있으며, 이 함수는 location 및 state 구조체를 HU2TransverseIsotropy에 전달합니다. 시뮬레이션 중에 FE 솔버는 location 및 state 구조체형 배열을 자동으로 계산하여 함수에 전달합니다.

ccoeffunc = @(location,state) HU2TransverseIsotropy(location,state,F);
specifyCoefficients(model,'m',0, ...
    'd',0, ...
    'c',ccoeffunc, ...
    'a',0, ...
    'f',[0;0;0]);

하중 적용 및 모델 해석하기

척추골의 축 방향 하중을 시뮬레이션하려면 아래쪽 표면을 고정하고 위쪽 표면에 하중을 아래 방향으로 가합니다. 분포 하중을 시뮬레이션하기 위해, 하중은 압력 형태로 적용됩니다.

경계 조건을 적용하기 위해 모델 형상의 면을 식별합니다. 이 예제에서는 pdegplot을 사용하여 위에서 만든 플롯의 시각적 검토를 통해 면이 식별되었습니다.

bottomSurfaceFace = 1;
topSurfaceFace = 250;

적용할 전체 힘을 뉴턴 단위로 지정합니다.

forceInput = -3000;

위쪽 표면의 면적을 추정합니다.

nf2=findNodes(model.Mesh,"region",Face=topSurfaceFace);
positions = model.Mesh.Nodes(:,nf2)';
surfaceShape = alphaShape(positions(:,1:2));
faceArea = area(surfaceShape);

압력 크기를 힘과 면적의 함수로 나타내어 계산합니다(단위: 파스칼).

inputPressure_Pa = forceInput/faceArea;

경계 조건을 적용합니다.

applyBoundaryCondition(model,"dirichlet",Face=bottomSurfaceFace,u=[0,0,0]);            
applyBoundaryCondition(model,"neumann",Face=topSurfaceFace,g=[0;0;inputPressure_Pa]);

모델을 해석합니다. 출력은 노드 변위와 그 기울기를 포함하는 StationaryResults (Partial Differential Equation Toolbox) 객체입니다.

Rs = solvepde(model);

결과 분석하기

모델 내의 축 방향 변위, 응력, 변형률을 플로팅하여 시뮬레이션 결과를 봅니다.

변위

pdeplot3D 함수를 사용하여 전체 모델의 축 방향 변위를 밀리미터 단위로 플로팅합니다.

Uz = Rs.NodalSolution(:,3)*10^3;
pdeplot3D(model,"ColorMapData",Uz)
clim([-0.15 0])
title("Axial Displacement (mm)")

Figure contains an axes object. The hidden axes object with title Axial Displacement (mm) contains 5 objects of type patch, quiver, text.

척추골 높이의 중간점에 위치한 횡단면 슬라이스의 변위를 플로팅합니다.

각 방향으로 1mm 간격으로 횡단면(xy) 슬라이스의 형상을 포함하는 사각 그리드를 만듭니다. x 및 y 벡터는 횡단면에서의 공간적 제한과 간격을 정의하고, z는 척추골 높이의 중간점을 제공하는 스칼라입니다. Xg, Yg 및 Zg 변수는 그리드 좌표를 정의합니다.

x = min(msh.Nodes(1,:)):0.001:max(msh.Nodes(1,:));
y = min(msh.Nodes(2,:)):0.001:max(msh.Nodes(2,:));
z = min(msh.Nodes(3,:))+0.5*(max(msh.Nodes(3,:))-min(msh.Nodes(3,:)));

[Xg,Yg] = meshgrid(x,y);
Zg = z*ones(size(Xg));

interpolateSolution (Partial Differential Equation Toolbox) 객체 함수를 사용하여 수직 축 방향 변위를 그리드 좌표에 보간합니다. 변위를 미터에서 밀리미터로 변환하고, 그리드 크기에 맞게 변위 벡터의 형태를 변경합니다.

U = interpolateSolution(Rs,Xg,Yg,Zg,3);
U = U*10^3;
Ug = reshape(U,size(Xg));

횡단면 슬라이스의 축 방향 변위를 플로팅합니다.

surf(Xg,Yg,Ug,LineStyle="none")
axis equal
grid off
view(0,90)
colormap("turbo")
colorbar
clim([-0.15 0])
title("Axial Displacement (mm)")

Figure contains an axes object. The axes object with title Axial Displacement (mm) contains an object of type surface.

응력

구조 해석에서 응력은 c 계수와 변위 기울기의 곱입니다. evaluateCGradient (Partial Differential Equation Toolbox) 객체 함수를 사용하여 횡단면 슬라이스의 그리드 좌표에서 수직 응력을 계산합니다.

[cgradx,cgrady,cgradz] = evaluateCGradient(Rs,Xg,Yg,Zg,3);

수직 축 방향 응력을 메가파스칼 단위로 변환하고, 그리드 크기에 맞게 응력 벡터의 형태를 변경합니다.

cgradz = cgradz*10^-6;
cgradzg = reshape(cgradz,size(Xg));

횡단면 슬라이스의 수직 축 방향 응력을 플로팅합니다.

surf(Xg,Yg,cgradzg,LineStyle="none");
axis equal
grid off
view(0,90)
colormap("turbo")
colorbar
clim([-10 1])
title("Axial Stress (MPa)")

Figure contains an axes object. The axes object with title Axial Stress (MPa) contains an object of type surface.

변형률

구조 해석에서 변형률은 변위 결과의 기울기입니다. evaluateGradient (Partial Differential Equation Toolbox) 객체 함수를 사용하여 횡단면 슬라이스의 그리드 좌표에서 수직 변형률을 추출합니다.

[gradx,grady,gradz] = evaluateGradient(Rs,Xg,Yg,Zg,3);

그리드 크기에 맞게 축 방향 변형률 벡터의 형태를 변경합니다.

gradzg = reshape(gradz,size(Xg));

횡단면 슬라이스의 축 방향 변형률을 플로팅합니다.

surf(Xg,Yg,gradzg,LineStyle="none");
axis equal
grid off
view(0,90)
colormap("turbo")
colorbar
clim([-0.008 0])
title("Axial Strain")

Figure contains an axes object. The axes object with title Axial Strain contains an object of type surface.

헬퍼 함수

HU2TransverseIsotropy 헬퍼 함수는 다음 단계를 사용하여 횡방향 등방성 물질 속성을 지정합니다.

scatteredInterpolant 객체 F를 사용하여 노드 위치에서 복셀 강도를 하운스필드 단위(HU)로 매핑합니다.
[1]의 선형 보정 방정식을 사용하여 HU를 CT 밀도로 변환합니다.
[2]의 밀도-탄성 관계식을 사용하여 CT 밀도를 축 방향의 탄성 계수 E3, 횡단면의 탄성 계수 E12, 체적 탄성 계수 G로 변환합니다. 축 방향의 푸아송 비 v3과 횡단면의 푸아송 비 v12 역시 [2]를 기준으로 할당됩니다.
elasticityTransOrthoC3D 헬퍼 함수를 사용하여 횡방향 등방성에 대한 6×6 컴플라이언스 행렬을 정의합니다.
SixMat2NineMat 및 SquareMat2CCoeffVec 헬퍼 함수를 사용하여 6×6 컴플라이언스 행렬, 즉 c 행렬을 Partial Differential Equation Toolbox에서 요구하는 벡터 형태로 변환합니다.

function ccoeff = HU2TransverseIsotropy(location,state,F)
    HU = F(location.x,location.y,location.z);
    rho = 5.2+0.8*HU;
    E3 = -34.7 + 3.230.*rho;
    E12 = 0.333*E3;
    v12 = 0.104;
    v3 = 0.381;
    G = 0.121*E3;

    ccoeff = zeros(45,length(location.x));
    for i = 1:length(location.x)
        cMatrix = elasticityTransOrthoC3D(E12(i),E3(i),v12,v3,G(i));
        nineMat = SixMat2NineMat(cMatrix);
        ccoeff(:,i) = SquareMat2CCoeffVec(nineMat);
    end
end

elasticityTransOrthoC3D 헬퍼 함수는 탄성 계수 E12와 E3, 체적 탄성 계수 G, 푸아송 비 v12와 v3을 기반으로 횡방향 등방성 선형 탄성 물질에 대한 6×6 컴플라이언스 행렬을 정의합니다. 헬퍼 함수는 컴플라이언스 행렬을 계산하기 전에 모든 탄성 계수 값을 메가파스칼 단위에서 파스칼 단위의 값으로 변환합니다.

function C = elasticityTransOrthoC3D(E12,E3,v12,v3,G)
    E12 = E12*10^6;
    E3 = E3*10^6;
    G = G*10^6;
    v_zp = (E3*v3)/E12;
    
    C = zeros(6,6);
    C(1,1) = 1/E12;
    C(1,2) = -v12/E12;
    C(1,3)= -v_zp/E3;
    C(2,1) = C(1,2);
    C(2,2) = 1/E12;
    C(2,3) = -v_zp/E3;
    C(3,1) = -v3/E12;
    C(3,2) = C(2,3);
    C(3,3) = 1/E3;
    C(4,4) = 1/(2*G);
    C(5,5) = 1/(2*G);
    C(6,6) = (1+v12)/E12;
    
    C=inv(C);

end

SixMat2NineMat 헬퍼 함수는 Voigt 표기법의 6×6 c 계수 행렬을 확장된 형태에 해당하는 9×9 행렬로 변환합니다.

function nineMat = SixMat2NineMat(sixMat)

    for i = 1:6
        nineVecs(i,:) = sixMat(i,[1 6 5 6 2 4 5 4 3]);
    end
    
    nineMat = [ nineVecs(1,:); ...
        nineVecs(6,:); ...
        nineVecs(5,:); ...
    
        nineVecs(6,:); ...
        nineVecs(2,:); ...
        nineVecs(4,:); ...
    
        nineVecs(5,:); ...
        nineVecs(4,:); ...
        nineVecs(3,:)];

end

SquareMat2CCoeffVec는 9×9 c 계수 행렬을 Partial Differential Equation Toolbox에서 요구하는 벡터 형태로 변환합니다. 이 헬퍼 함수는 c Coefficient for specifyCoefficients (Partial Differential Equation Toolbox)의 "3N(3N+1)/2개 요소를 가진 열 벡터 c, 3차원 시스템" 경우에 해당하는, 요소를 45개 가진 벡터를 생성합니다.

function c45Vec = SquareMat2CCoeffVec(nineMat)

    C11 = [nineMat(1,1) nineMat(1,2) nineMat(2,2) nineMat(1,3) nineMat(2,3) nineMat(3,3)];
    C12 = [nineMat(1,4) nineMat(2,4) nineMat(3,4) nineMat(1,5) nineMat(2,5) nineMat(3,5) ...
        nineMat(1,6) nineMat(2,6) nineMat(3,6)];
    C13 = [nineMat(1,7) nineMat(2,7) nineMat(3,7) nineMat(1,8) nineMat(2,8) nineMat(3,8) ...
        nineMat(1,9) nineMat(2,9) nineMat(3,9)];
    C22 = [nineMat(4,4) nineMat(4,5) nineMat(5,5) nineMat(4,6) nineMat(5,6) nineMat(6,6)];
    C23 = [nineMat(4,7) nineMat(5,7) nineMat(6,7) nineMat(4,8) nineMat(5,8) nineMat(6,8) ...
        nineMat(4,9) nineMat(5,9) nineMat(6,9)];
    C33 = [nineMat(7,7) nineMat(7,8) nineMat(8,8) nineMat(7,9) nineMat(8,9) nineMat(9,9)];
    c45Vec  = [C11 C12 C22 C13 C23 C33]';

end

참고 문헌

[1] Turbucz, Mate, Agoston Jakab Pokorni, György Szőke, Zoltan Hoffer, Rita Maria Kiss, Aron Lazary, and Peter Endre Eltes. “Development and Validation of Two Intact Lumbar Spine Finite Element Models for In Silico Investigations: Comparison of the Bone Modelling Approaches.” Applied Sciences 12, no. 20 (January 2022): 10256. https://doi.org/10.3390/app122010256.

[2] Unnikrishnan, Ginu U., and Elise F. Morgan. “A New Material Mapping Procedure for Quantitative Computed Tomography-Based, Continuum Finite Element Analyses of the Vertebra.” Journal of Biomechanical Engineering 133, no. 7 (July 1, 2011): 071001. https://doi.org/10.1115/1.4004190.

참고 항목

도움말 항목

의료 영상 좌표계
Label 3-D Medical Image Using Medical Image Labeler
Create STL Surface Model of Femur Bone for 3-D Printing
c Coefficient for specifyCoefficients (Partial Differential Equation Toolbox)