Changeset 86


Ignore:
Timestamp:
Sep 4, 2009, 2:50:44 PM (10 years ago)
Author:
tjniskan
Message:
  • Fixed a image scaling bug in fat analysis -plugin

M plugins/fat_analysis.m
M aedes_revision.m

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • aedes_revision.m

    r85 r86  
    9393% bash-script every time it is called so that this file "aedes_revision.m" is
    9494% always in the list of committed files. DO NOT EDIT THE NEXT LINE!!!
    95 % - SVN Hook -
     95% - Svn Hook -
  • plugins/fat_analysis.m

    r83 r86  
    108108end
    109109
     110% A brute force -method for correct DC artifacts, may not always work
     111data_water = DATA{1}.FTDATA(:,:,:,2);
     112data_fat = DATA{1}.FTDATA(:,:,:,1);
     113[mx,DC_Ind_fat]=max(max(data_fat(:,:,data_sz(3)/2),[],2));
     114[mx,DC_Ind_water]=max(max(data_water(:,:,data_sz(3)/2),[],2));
     115
     116% Check that the found noise peak is near the center
     117if DC_Ind_fat < data_sz(1)/2-5 || DC_Ind_fat > data_sz(1)/2+5
     118  warning('Largest noise peak not near the center line in fat data! Check images and results.');
     119  DC_Ind_fat=data_sz(1)/2+1;
     120end
     121if DC_Ind_water < data_sz(1)/2-5 || DC_Ind_water > data_sz(1)/2+5
     122  warning('Largest noise peak not near the center line in water data! Check images and results.');
     123  DC_Ind_water=data_sz(1)/2+1;
     124end
     125data_fat(DC_Ind_fat,:,:) = mean(data_fat([DC_Ind_fat-1 DC_Ind_fat+1],:,:));
     126data_water(DC_Ind_water,:,:) = mean(data_water([DC_Ind_water-1 DC_Ind_water+1],:,:));
     127
     128
    110129%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    111130%% FAT analysis -----------------------------------------
     
    113132
    114133%% Set the default noise ROI
    115 noise_fat = 2*squeeze(mean(mean(DATA{1}.FTDATA(data_sz(1)-10:data_sz(1),1:data_sz(2),:,1),1),2));
    116 noise_water = 2*squeeze(mean(mean(DATA{1}.FTDATA(data_sz(1)-10:data_sz(1),1:data_sz(2),:,2),1),2));
    117 
    118 mean_fat = squeeze(mean(mean(DATA{1}.FTDATA(:,:,:,1),1),2));
    119 mean_water = squeeze(mean(mean(DATA{1}.FTDATA(:,:,:,2),1),2));
     134% noise_fat = 2*squeeze(mean(mean(DATA{1}.FTDATA(data_sz(1)-10:data_sz(1),1:data_sz(2),:,1),1),2));
     135% noise_water = 2*squeeze(mean(mean(DATA{1}.FTDATA(data_sz(1)-10:data_sz(1),1:data_sz(2),:,2),1),2));
     136noise_fat = 2*squeeze(mean(mean(data_fat(1:30,1:40,:),1),2));
     137noise_water = 2*squeeze(mean(mean(data_water(1:30,1:40,:),1),2));
     138
     139mean_fat = squeeze(mean(mean(data_fat(:,:,:),1),2));
     140mean_water = squeeze(mean(mean(data_water(:,:,:),1),2));
    120141
    121142fat_th = noise_fat+mean_fat;
     
    124145water_th = water_th(:);
    125146
    126 %z_range = 20:78;
    127 %This value comes from aedes: see which slices correspond to right
    128 %anatomical area and have tolerable level of noise/artefacts
    129 
    130 % sum_fat = 0;
    131 % sum_h2o = 0;
    132 % for ind=z_range
    133 %     sum_fat=sum_fat + sum(sum((DATA{1}.FTDATA(:,:,ind,1)>fat_th(ind))));
    134 %     sum_h2o=sum_h2o + sum(sum((DATA{1}.FTDATA(:,:,ind,2)>water_th(ind))));
    135 % end
    136 
    137147% Construct matrices for threshold comparison
    138148fat_ind = repmat(permute(fat_th(z_range),[3 2 1]),...
    139   [size(DATA{1}.FTDATA,1) size(DATA{1}.FTDATA,2) 1]);
     149  [data_sz(1) data_sz(2) 1]);
    140150water_ind = repmat(permute(water_th(z_range),[3 2 1]),...
    141   [size(DATA{1}.FTDATA,1) size(DATA{1}.FTDATA,2) 1]);
     151  [data_sz(1) data_sz(2) 1]);
    142152
    143153% Calculate number of water and fat voxels
    144 nVoxFat = numel(find(DATA{1}.FTDATA(:,:,z_range,1)>fat_ind));
    145 nVoxWater = numel(find(DATA{1}.FTDATA(:,:,z_range,2)>water_ind));
     154nVoxFat = numel(find(data_fat(:,:,z_range)>fat_ind));
     155nVoxWater = numel(find(data_water(:,:,z_range)>water_ind));
    146156
    147157% Calculate fat and water concentrations
     
    214224
    215225if ispc
    216   set(h_text,'fontsize',10);
     226  set(h_text,'fontsize',8);
    217227else
    218228  set(h_text,'fontsize',8);
     
    274284end
    275285
    276 % Correct DC artifact
    277 data_water = DATA{1}.FTDATA(:,:,:,2);
    278 data_fat_shifted = DATA{1}.FTDATA(:,:,:,1);
    279 [mx,DC_Ind_fat]=max(max(DATA{1}.FTDATA(:,:,size(DATA{1}.FTDATA,3)/2,1)));
    280 [mx,DC_Ind_water]=max(max(DATA{1}.FTDATA(:,:,size(DATA{1}.FTDATA,3)/2,2)));
    281 data_fat_shifted(DC_Ind_fat,:,:) = mean(data_fat_shifted([DC_Ind_fat-1 DC_Ind_fat+1],:,:));
    282 data_water(DC_Ind_water,:,:) = mean(data_water([DC_Ind_water-1 DC_Ind_water+1],:,:));
     286
     287% A brute force -method for correct DC artifacts, may not always work
     288% data_water = DATA{1}.FTDATA(:,:,:,2);
     289% data_fat_shifted = DATA{1}.FTDATA(:,:,:,1);
     290% [mx,DC_Ind_fat]=max(max(DATA{1}.FTDATA(:,:,size(DATA{1}.FTDATA,3)/2,1),[],2));
     291% [mx,DC_Ind_water]=max(max(DATA{1}.FTDATA(:,:,size(DATA{1}.FTDATA,3)/2,2),[],2));
     292% data_fat_shifted(DC_Ind_fat,:,:) = mean(data_fat_shifted([DC_Ind_fat-1 DC_Ind_fat+1],:,:));
     293% data_water(DC_Ind_water,:,:) = mean(data_water([DC_Ind_water-1 DC_Ind_water+1],:,:));
    283294
    284295%% Calculate Fat shift ------------------
     
    288299fat_shift = 680/sw;
    289300
    290 data_fat_shifted = interp1(1:size(data_fat_shifted,1),...
    291   data_fat_shifted,...
    292   (1:size(data_fat_shifted,1))+fat_shift,'spline',0);
     301data_fat = interp1(1:size(data_fat,1),...
     302  data_fat,...
     303  (1:size(data_fat,1))+fat_shift,'spline',0);
     304
     305% If there is a huge DC artifact in the data the spile interpolation tends
     306% to go crazy near the artifacts and produce huge negative values...
     307data_fat(data_fat<0) = 0;
     308
     309% Close waitbar
    293310close(h)
    294311% ---------------------------------------
    295312
    296313% Try to set the intensity levels to reasonable values
    297 tmp_water=data_water((size(data_water,1)/2-20):(size(data_water,1)/2+20),...
    298   (size(data_water,2)/2-20):(size(data_water,2)/2+20),:);
    299 tmp_fat=data_fat_shifted((size(data_fat_shifted,1)/2-20):(size(data_fat_shifted,1)/2+20),...
    300   (size(data_fat_shifted,2)/2-20):(size(data_fat_shifted,2)/2+20),:);
    301 data_fat_shifted = data_fat_shifted./max(tmp_fat(:));
     314%tmp_water=data_water;
     315tmp_water = data_water((size(data_water,1)/2-30):(size(data_water,1)/2+30),...
     316  (size(data_water,2)/2-30):(size(data_water,2)/2+30),:);
     317tmp_fat = data_fat((size(data_fat,1)/2-30):(size(data_fat,1)/2+30),...
     318  (size(data_fat,2)/2-30):(size(data_fat,2)/2+30),:);
     319data_fat = data_fat./max(tmp_fat(:));
    302320data_water = data_water./max(tmp_water(:));
    303 data_fat_shifted(data_fat_shifted>1) = 1;
     321data_fat(data_fat>1) = 1;
    304322data_water(data_water>1) = 1;
    305323
    306324% Construct an RGB image of water and fat content
    307 RGB_im = cat(4,data_fat_shifted,...
    308   zeros(size(data_fat_shifted)),...
     325RGB_im = cat(4,data_fat,...
     326  zeros(size(data_fat)),...
    309327  data_water);
    310328
     
    328346end
    329347close(h)
    330 
    331 %{
    332 
    333 %%%%%%%%%%%%%% FAT shift %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    334 
    335 sw =  65627.5635767/siz(1) ; % per pixel 65627.5635767 FOV 40x40, 65789.4736842 FOV 50x50
    336 shift = 680/sw %in pixels
    337 
    338 im_fat2 = zeros(siz(1),siz(2),siz(3));
    339 
    340 for ind2=1:siz(3)
    341     for ind=1:siz(2)
    342             im_fat2(:,ind,ind2) =  interp1(1:siz(1),im_fat(:,ind,ind2),(1:siz(1))+shift,'spline',0);
    343     end
    344 end
    345 
    346 
    347 %%%%%%%%%%%%%%%%%%%% RGB image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    348 
    349 %RGB_im=cat(4,im_fat,im_h2o, zeros(siz(1),siz(2),siz(3)));
    350 
    351 mask_h2o=ones(siz(1),siz(2));
    352 mask_fat=ones(siz(1),siz(2));
    353 r1_h2o = [(siz(1)/2-2):(siz(1)/2+2)]; %to exclude the middle (DC offset)
    354 r1_fat = [(siz(1)/2-5):(siz(1)/2-1)];
    355 r2 = [(siz(2)/2-2):(siz(2)/2+2)];
    356 mask_h2o(r1_h2o,r2) = 0.001;
    357 mask_fat(r1_fat,r2) = 0.001;
    358 
    359 im_fat3=im_fat;
    360 im_h2o3=im_h2o;
    361 
    362 for ind=1:siz(3)
    363     im_fat3(:,:,ind)=im_fat2(:,:,ind)/(max(max(mask_fat.*squeeze(im_fat2(:,:,ind)))));
    364     im_h2o3(:,:,ind)=im_h2o(:,:,ind)/(max(max(mask_h2o.*squeeze(im_h2o(:,:,ind)))));
    365 end
    366 
    367 %   im_fat3(:,:,63)=im_fat2(:,:,65)/(1*max2(mask_fat.*squeeze(im_fat2(:,:,65))));
    368 %   im_h2o3(:,:,63)=im_h2o(:,:,65)/(1*max2(mask_h2o.*squeeze(im_h2o(:,:,65))));
    369 
    370 im_fat3(im_fat3>1)=1;
    371 im_h2o3(im_h2o3>1)=1;
    372 im_fat3(im_fat3<0)=0;
    373 im_h2o3(im_h2o3<0)=0;
    374 
    375 RGB_im=abs(cat(4,im_fat3, zeros(siz(1),siz(2),siz(3)),im_h2o3));
    376 imagesc([[squeeze(abs(RGB_im(:,:,40,:))),squeeze(abs(RGB_im(:,:,50,:))),squeeze(abs(RGB_im(:,:,60,:)))];[squeeze(abs(RGB_im(:,:,70,:))),squeeze(abs(RGB_im(:,:,80,:))),squeeze(abs(RGB_im(:,:,90,:)))]])
    377 axis image
    378 
    379 clear im_fat
    380 clear im_h2o
    381 clear im_fat2
    382 
    383 %thresh_scale = 5;
    384 %fatnoise_rel=squeeze(mean(mean(im_fat3(noiseroix, noiseroiy,:),1),2));
    385 %h2onoise_rel=squeeze(mean(mean(im_h2o3(noiseroix, noiseroiy,:),1),2));
    386 
    387 %fat_rel = 0;
    388 %h2o_rel = 0;
    389 %for ind=z_range
    390 %    fat_rel = fat_rel + sum(sum(sum(im_fat3(im_fat3(:,:,ind)>(thresh_scale*fatnoise_rel(ind))))))* density_fat/1000; %- siz(1)*siz(2)*siz(3)*mean(fatnoise_rel)/4*density_fat/1000
    391 %    h2o_rel = h2o_rel + sum(sum(sum(im_h2o3(im_h2o3(:,:,ind)>(thresh_scale*h2onoise_rel(ind))))))* density_h2o/1000; %- siz(1)*siz(2)*siz(3)*mean(h2onoise_rel)/4*density_h2o/1000
    392 %end
    393 
    394 %fat_rel
    395 %h2o_rel
    396 %fat_percent_rel = fat_rel/(fat_rel+h2o_rel)
    397 
    398 for ind=1:siz(3)
    399     imwrite(squeeze(abs(RGB_im(:,:,ind,:))),['P:\nmr\NEW_Johanna\data\mousefat\Fatmouse_m09_', num2str(ind),'.tif'],'tif');
    400 end
    401 
    402 %}
    403 
    404 
    405 
    406 
Note: See TracChangeset for help on using the changeset viewer.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008