Tuesday, 21 July 2015

Wrote the following code to calculate the gradient of an isosurface of tumour to give a quantitative measure of shape:


%% gradient of a slice of tumour
% if not already defined, make Tumour from cstore
% for i = 1:201
% Tumour(:,:,:,i) =cstore{i};
% end
% then take time point t
% Tumour(t) = Tumour(:,:,:,t)
t = 200;
Slice = Tumour(:,:,40,t); %slice at z = 40, t = 125

%Manually make isosurface matrix
%RR = isosurface(Tumour80,2e5); then we have sructure RR.V (1613x3 double)
%and RR.F
RR = zeros(102,102);
ISO = Slice >2e5;
RR = RR + ISO;
%calculate gradient
%%Anisotropy Measure using gradient integral
%A quantitative measure of tumour anisotropy,
%1/(2*pi*R)*Integraldxdy* absolute value of grad phi
%discretised version:
%Sigma(ijk) abs(grad c(r)) =
%sqrt((partialc/partialx)^2+(partialc/partialy)^2+(partialc/partialz))dV
%partialc/partialx = ((Ci+1,j,k - Ci-1,j,k)/2deltax)

%practice code

dx = 1;
for i = 1:99, j = 1:100;
    dcdxsq =(RR(i+1,j) - RR(i,j)).^2;
end
for i=1:100, j = 1:99;
    dcdysq=(RR(i,j+1) - RR(i,j)).^2;
end
GridSum = sqrt(sum(sum(dcdxsq)) + sum(sum(dcdysq)));
figure
imagesc(RR)

GG=gradient(RR);
Aniso2 = sqrt(sum(sum(GG.^2)))
str = sprintf('Tumour surface = %f',Aniso2);
title(str);
%ThreeD version

TSurf = Tumour(:,:,:,t); %tumour matrix, t = 130
TIso = TSurf >2e5;
TTT = zeros(102,102,62);
TIso = TIso + TTT;
figure
fv = isosurface(TSurf,2e5);
p1 = patch(fv,'FaceColor','red','EdgeColor','none');
view(3)
daspect([1,1,1])
axis([0 100 0 100 0 60])
camlight
camlight(-80,-10)
lighting gouraud
title('Triangle Normals')
IsoGradient = gradient(fv.vertices);
AbsGrad3 = sqrt(sum(sum(IsoGradient.^2)))
str = sprintf('Tumour isosurface with Anisotropy = %f', AbsGrad3);
title(str)
%total tumour cells inside isosurface
ZZ = find(TSurf>2e5);
ZZZ= sum(sum(sum(ZZ)));

stf = sprintf('total tumour cells inside isosurface = %f', ZZZ)
str
Aniso2
AbsGrad3

Using this, I tabulated and plotted the data from a range of anisotropies using this code
%graphing anisotropy\
AnisoData = ...
[  NaN,      0,       1,            10,       15,         30,          45,       60,    100; ...
   NaN,    NaN,     NaN,           NaN,      NaN,         NaN,        NaN,      NaN,    NaN; ...
    40,      0,       0,             0,        0,          0,           0,        0,      0; ...
    50,      0,       0,          50.8308,  50.8358,     50.8358,    50.8572,   50.8621,   50.8690; ...
    60,   531.8208, 273.9941,    272.4435,  274.2529,   274.2529,   272.1032,   272.1409,  272.5392; ...
    70,   978.4913, 498.2811,    511.0714,  511.8223,   511.8223,   514.6845,   515.4931,  521.0318; ...
    80, 1.3849e+03, 764.0871,    784.7974,  787.4831,   787.4831,   802.9823,   802.8959,  809.3151; ...
    90, 1.8490e+03, 1.0832e+03, 1.1284e+03, 1.1359e+03, 1.1359e+03, 1.1695e+03, 1.1801e+03, 1.1909e+03; ...
   100, 2.3453e+03, 1.4825e+03, 1.5795e+03, 1.5993e+03, 1.5993e+03, 1.6569e+03, 1.6642e+03,1.6791e+03; ...
   110, 2.9061e+03, 1.9680e+03, 2.1226e+03, 2.1498e+03, 2.1498e+03, 2.2122e+03, 2.2211e+03,2.2375e+03; ...
   120, 3.4490e+03, 2.4469e+03, 2.6406e+03, 2.6727e+03, 2.6727e+03, 2.7469e+03, 2.7582e+03, 2.7812e+03; ...
   130, 3.9770e+03, 2.9079e+03, 3.1535e+03, 3.2000e+03, 3.2000e+03, 3.2964e+03, 3.3147e+03, 3.3333e+03; ...
   140, 4.4230e+03, 3.3718e+03, 3.6241e+03, 3.6663e+03, 3.6663e+03, 3.7589e+03, 3.7653e+03, 3.7858e+03; ...
   150, 4.6209e+03, 3.8110e+03, 4.0923e+03, 4.1365e+03, 4.1365e+03, 4.2477e+03, 4.2624e+03, 4.2996e+03; ...
   160, 4.8136e+03, 4.2563e+03, 4.5334e+03, 4.5946e+03, 4.5946e+03, 4.6870e+03, 4.6991e+03, 4.7366e+03; ...
   170, 4.9855e+03, 4.6604e+03, 4.9125e+03, 4.9501e+03, 4.9501e+03, 5.0425e+03, 5.0592e+03, 5.0788e+03; ...
   180, 5.1144e+03, 5.0263e+03, 5.1811e+03, 5.2162e+03, 5.2162e+03, 5.2552e+03, 5.2698e+03, 5.2814e+03; ...
   190, 5.2438e+03, 5.2282e+03, 5.2929e+03, 5.2902e+03, 5.2902e+03, 5.2949e+03, 5.2953e+03, 5.3019e+03; ...
   200, 5.2867e+03, 5.2929e+03, 5.2865e+03, 5.2865e+03, 5.2865e+03, 5.2865e+03, 5.2865e+03, 5.2865e+03];
X = AnisoData(3:end,1)*13.82;
Y0 = AnisoData(3:end,2);
Y1 = AnisoData(3:end,3);
Y10 = AnisoData(3:end,4);
Y15 = AnisoData(3:end,5);
Y30 = AnisoData(3:end,6);
Y45 = AnisoData(3:end,7);
Y60 = AnisoData(3:end,8);
Y100 = AnisoData(3:end,9);
plot(X,Y0, X,Y1, X,Y10, X,Y15, X,Y30, X,Y45, X,Y60, X,Y100)
axis([0 3000 0 5.5e3])
legend('y = rf0','y = rf1','y = rf10','y = rf15','y = rf30','y = rf45','y = rf60','y = rf100')
title('Plot of gradient of isosurface for varying anisotropy')
xlabel('time of tumour growth (days)')
ylabel('Abs(grad(isosurface))')


and produced some preliminaty graphs, which show the data behaves as expected and this may be a good measure of anisotropy of tumour surface.

Tomorrow I will explore the sum of tumour cell concentration within an isosurface as the anisotropy is varied. 

Then on to modelling resection.

Meeting with Z yesterday, plan is to explore virtual surgery this year, talk about experiments with mouse models next year and fully get to grips with what we have here and get my skills up a bit before moving on to more realistic simulations involving tissue mechanics. 
As Zoltan said, my skills are probably not up to tackling this kind of problem just yet, so this sems like a very good plan to me.


No comments:

Post a Comment