% dtry37 Make possible to choose 64, 128 or 256 when rect fov. % NOTE: it complains on rect 128 case, but does it right. % diff1 Modify for diffusion sequences; no distcor for now. % Change matrix size, take out loopcor. % diff2 Put in /working directory in case alliance goes down. % diff3 Take out the ghost correction that uses calibration. % diff4 Put in Shy's ghost correction routine. % diff5 Add Haselgrove distortion correction. % Threshold of 20%, works better than 15%. % diff6 Was correcting mag and trans in WRONG direction. 84 is phase dir. % diff8 Ghost correct raw images, then average them before dist. correct. % diff9 Add a median filter before distortion correct % diff10 Fix slope and start, as in difftune7.m % diff11 Use steepest gradient as edge. % Modify so translation can't jump big. Plan to use on phantoms. % diff14 Modify so don't need so much memory; Fix one frame at a time. % Takes 10 min till ready, as opposed to 17 min for diff11. % diff16 Fit the magnification to a line to take out noise. % diff17 Scale images AFTER the ft, and write individually. % bip Make a version for bipolar data, needs no correction. % bip2 Save the reference image set also; % Redo for Robert, Ref, Anisotropy, BY. % bip3 Read in displace matrix, and do distortion correction. % Matrices are different sizes. Do correcting on 256 matrix. % FIX FOV ASPECT RATIO IN THE DIFFUSION IMAGES. WAS 6% OFF!! % Use bilinear interp. in imresize; 'nearest neighbor' led % to pixelated look that caused ridges in ramp image. % WORKING VERSION FOR AXIAL DIFFUSION. DOES DISTOR CORRECTION patfile = '/working/fmri/bruce2/kspace.20'; newname= '/working/fmri/bruce2/difflow'; disname= '/working/fmri/bruce2/displace'; slices=10; % The number of slices in ghost calib and patfile averages=10; % The number of averages in patfile frames=8; % Code is hardwired for 8 diffusion frames. % Sequence does a dummy to set TR, then Bo, then 6 grads. % Normally firstl= 2 lastsl=11 for low set and % firstsl=12 lastsl=21 for high one. firstsl=2; % First and last slices of displace to use in distor. corr. lastsl=11; pctlevel = 0.12; b_fact = 1070; % Calibrated this value Experimentally. See page 31 book 3. fermifilter= 'y'; size=256; % NOTE: assume size is 256; and use 272 for aspect ratio of % 272/256 = 1.0625. Can make general for 128 later. fermi1 = 0.8; % Constants to control fermi filter fermi2 = 50; fermi = zeros(90,84); f3 = zeros(90,84); test2=zeros(90,84); % Set up the fermi filter matrix % NOTE: This uses half the image dimensions. for j=1:45 for k=1:42 rj = j/45; rk = k/42; r = rj*rj + rk*rk; filt = 1.0 / ( 1.0 + exp((r-fermi1)*fermi2)); fermi(45+j,42+k) = filt; fermi(45-j+1,42+k) = filt; fermi(45+j,42-k+1) = filt; fermi(45-j+1,42-k+1) = filt; end end % First long word is length of header % Second long word is length of data in bytes % Divide by four to get number of complex pixels % Use b1 array with int16 read to get high slice numbers % PATIENT DATA ______________________________________________________ fid=fopen(patfile,'r','ieee-be'); a=fread(fid,2,'long'); b1=zeros(slices,averages,a(1)/2); F1=zeros(slices,averages,a(2)/4); FF=zeros(slices,averages,a(2)/8); sumavg = zeros(90,84,slices,frames); % Do a dummy read for the first frame, set of slices, to throw out. % Use it to choose boundaries for ghosting test2 = zeros(90,84); for i=1:slices b1(i,1,:)=fread(fid,a(1)/2,'uint16').'; F1(i,1,:)=fread(fid,a(2)/4,'float32').'; FF(i,1,:)=F1(i,1,1:2:end)+1i*F1(i,1,2:2:end); test2 = test2 + reshape(FF(i,1,:),90,84); end % Pick boundaries for ghosting f3 = test2.'; [s1,s2]=size(f3); odd=zeros(s1,s2);; odd(1:2:s1,:)=f3(1:2:s1,:); new_odd=fftshift(ifft2(odd)); imagesc(abs(new_odd)); disp('Click on upper left and lower right no-ghost corners'); [x,y]=ginput(2); % Start a loop to read patient data for frame=1:(frames-1) avg = zeros(90,84,slices,averages); disp(sprintf('Reading, ghost correcting, averaging frame %d',frame)); for j = 1:averages for i = 1:slices b1(i,j,:)=fread(fid,a(1)/2,'uint16').'; F1(i,j,:)=fread(fid,a(2)/4,'float32').'; FF(i,j,:)=F1(i,j,1:2:end)+1i*F1(i,j,2:2:end); fr = floor ( (b1(i,j,12)/slices) - 0.01) + 1; sl = b1(i,j,12) - ( (fr-1)*slices); avg(:,:,sl, j)=reshape(FF(i,j,:),90,84); %Fermi filter each image if (fermifilter == 'y') avg(:,:,sl,j)=avg(:,:,sl,j).* fermi;end % Ghost correct the patient data f3=avg(:,:,sl,j).'; [s1,s2]=size(f3); even=zeros(s1,s2); odd=even; odd(1:2:s1,:)=f3(1:2:s1,:); even(2:2:s1,:)=(f3(2:2:s1,:)); new_odd=fftshift(ifft2(odd)); new_even=fftshift(ifft2(even)); theta=angle(new_even)-angle(new_odd); theta(1:floor(min(y))-1,:)=0; theta((ceil(max(y))+1):end,:)=0; theta(:,1:floor(min(x)-1))=0; theta(:,(ceil(max(x))+1):end)=0; theta(abs(new_odd)pi)=theta(theta>pi)-2*pi; theta(theta<-pi)=theta(theta<-pi)+2*pi; x_ky=ifft(odd+even,[],2); theta_x=zeros(1,s2); for k=1:s2 data_x(k)=sum(theta(:,k)~=0); if data_x(k)>10 theta_x(k)=sum(theta(:,k).*(theta(:,k)~=0))/((data_x(k))*2); end end for k=1:s1 if rem(k,2)==1 x_ky(k,:)=x_ky(k,:).*fftshift(exp(1i*theta_x)); else x_ky(k,:)=x_ky(k,:).*fftshift(exp(-1i*theta_x)); end end avg(:,:,sl,j)=fft(x_ky,[],2).'; % % FT the avg array avg(:,:,sl,j) = abs(fftshift(ifft2(avg(:,:,sl,j)))); end % End of slices end % End of averages % AVERAGE THE PATIENT IMAGES _____________________________________________ % Sum the corrected averages. % sumavg will have entries for 2-8; 1 is blank. for i=1:slices for j=1:averages sumavg(:,:,i,frame+1) = sumavg(:,:,i,frame+1) + avg(:,:,i,j)/averages; end end end % End of frame loop fclose(fid); clear b1; clear FF; clear F1; clear avg; % Change matrix to size x size ___________________________________________ bigavg=zeros(size,size,slices,frames); test4=zeros(272,256); test5=test4; begin1 = (272 - 90 ) / 2; begin2 = (256 - 84 ) / 2; ahalf=(272-256)/2; % Write fft of sumavg into center of larger array test4 for i=1:slices for j=2:frames test4(begin1+1:begin1+90,begin2+1:begin2+84) ... = fftshift(fft2(sumavg(:,:,i,j))); test5=ifft2(test4); bigavg(:,:,i,j) = abs(fliplr(flipud(test5((ahalf+1):(272-ahalf),:)))); end end clear sumavg; newavg=zeros(size,size,slices,frames); % % BRING IN DISPLACE AND RESIZE TO 256x256----------------------- disp('Starting to read displace matrix'); displace=zeros(132,68,22); bigdis=zeros(256,256,22); test3=zeros(554,256); aspect=554; final2=zeros(aspect,256); begin1=(aspect-132)/2; begin2=(size-132)/2; half = (aspect-size)/2; % Scale the matrix by 256/68 to shift more pixels, since new pixels are smaller. fid3=fopen(disname,'r','native'); for i=1:22 f1=fread(fid3, 8976, 'int16'); displace(:,:,i)=(256/68)*reshape(f1,132,68); end fclose(fid3); % Change size, then cut out middle. Do flips to match what done above. % Negative sign is necessary because of the different orientation. kk=1; for i=firstsl:lastsl test3 = imresize(displace(:,:,i),[554 256],'bilinear'); bigdis(:,:,kk) = -(fliplr(flipud(test3((half+1):(aspect-half),:)))); kk=kk+1; end % Do the ramp matrix. ramp=zeros(256,256,10); disp('Create ramp matrix'); for slice=1:slices for line=1:256 for j=2:255 [p,s]=polyfit([1,2,3],[bigdis(line,j-1,slice),bigdis(line,j,slice),bigdis(line,j+1,slice)],1); ramp(line,j,slice) = 1 + ( p(1)/10 ); end end end % Ramp and displace matrices are done. Start correcting data. % Size was 671. n-1*10 then add 1. 2550 +1 =2551 % In this code shifty was original, noshift was new % Need to create another huge matrix for the results for frame=2:frames disp(sprintf('Correcting frame %d',frame)); for slice=1:slices for line=1:256 newreal=repmat(0,1,2551); if (max(abs(bigdis(line,:,slice))) > 0) for j=2:255 % Divide ramp value by 10 since the bigdis displace matrix is scaled by 10 % NOTE: moved the operation up to ramp formation line % Use 5 since want half of total for index pixsize = round( 5 * ramp(line,j,slice)); if (pixsize > 0) pixsig = bigavg(line,j,slice,frame); elseif (pixsize <= 0) pixsig=0; end % Throw out points that go out of image address = (j-1)*10+1+floor(bigdis(line,j,slice)); if ( (address <=2551)& (address >= 1) ) % Spread out the pixel values in size for s=-(pixsize-1):pixsize if ( ((address+s) < 2551) & ((address+s) > 1) ) newreal(address+s) = newreal(address+s) + pixsig; end end % End of s loop end % End of if end % End of j % Sum up 10 points for each new pixel % WHERE 13 and 4 COME FROM? point 2 maps to 11; 20-13=7, and 20-4=16, so these % are the 10 points symmetrically about the 2nd pixel. for j=2:255 newavg(line,j,slice, frame) = sum(newreal( (10*j-13):(10*j-4) )); end else % If displace line is all zeros, just copy newavg(line,:,slice,frame) = 10*bigavg(line,:,slice, frame); end % End of if on displacement end % End of line loop end % End of slice loop end % End of frame loop % CALCULATE IMAGES ----------------------------------------------------------------- %clear bigavg; %clear bigdis; disp('Calculate Diffusion Images'); % Declare matrices for below testa=zeros(size,size); lnimg = zeros(size,size,slices,frames); d = lnimg; adc = zeros(size,size,slices); isowt2=adc; mask=adc; superimp=adc; aniso=adc; testb=adc; ref=adc; BY=adc; BX=adc; % Now compute all the computed images. % newavg has Bo in 2, 3-BY 4-BX 5-BZ single grads, % 6-8 double grads for each slice. % scale newavg so max is 2000. On 5 mm slices of me, max was 2.1 % Bigavg defined for frames 2-8. Frame 1 is the dummy. bigmax = max(max(max(max(newavg)))); newavg = 2000*newavg/bigmax; % Put ref images in matrix ref for i=1:slices ref(:,:,i)=newavg(:,:,i,2); end % Put BY in matrix BY for i=1:slices BY(:,:,i)=newavg(:,:,i,3); end % Put BX in matrix BY for i=1:slices BX(:,:,i)=newavg(:,:,i,4); end % Put 3rd % LNIMG; defined for frames 2 - 8, since newavg(1)=0; lnimg (newavg ~= 0) = log( 2 * newavg (newavg~= 0)); % D IMAGES; defined for frames 3-8; uses 2 as reference. for j=3:frames d(:,:,:,j) = 10000. * (lnimg(:,:,:,2)- lnimg(:,:,:,j) ) / b_fact; end clear lnimg; % ISOWT2 IMAGES for i=1:slices for j=3:5 isowt2(:,:,i) = isowt2(:,:,i) + newavg(:,:,i,j).*newavg(:,:,i,j); end end % Take off the scale by 10 that was on isowt2 in diffuse4. % The ranges of values are different here since I scale to 2000 later. isowt2 = sqrt(isowt2) ; % MASK : uses 2-5, since that is used in diffuse4.c . for j=2:5 superimp(:,:,:) = superimp(:,:,:) + newavg(:,:,:,j); end %clear newavg; % JNL. New software did not let me do this as 3D. %superimp = medfilt2(superimp, [21 21]); for i=1:slices superimp(:,:,i)=medfilt2(superimp(:,:,i), [21 21]); end supermax = max(max(max(superimp))); level = pctlevel * supermax; mask(superimp > level) = 1; % ADC ; one for each slice. Use values of d > 0 to keep out negatives. for j=3:5 testb = squeeze(d(:,:,:,j)); % remove last singleton dim. adc(testb > 0) = adc(testb > 0) + 0.3333333333 * testb(testb>0); end % Mask the outer edges of adc so the max routine below is not confused % by noise adc = mask .* adc; topd=repmat(0,1,slices); for i=1:slices topd(i) = max(max( adc(:,:,i) ) ); end % ANISTROPY for j=3:8 testb=squeeze(d(:,:,:,j)); aniso = aniso + (adc - testb) .* (adc - testb); end adc1 = adc; % Copy this to use for scaling below for i=1:slices testa=adc1(:,:,i); testa (testa < (0.3*topd(i))) = 0.3*topd(i); adc1(:,:,i) = testa; end aniso = mask .* sqrt (aniso) ./ adc1; % Now that you are done calculating, multiply adc by mask adc = adc .* mask; % Scale computed images images so the max is 2000 % scale using filtered versions since otherwise noise points can end up being % the max in the volume for i=1:slices testb(:,:,i)=medfilt2(isowt2(:,:,i), [21 21]); end bigmax = max(max(max(testb))); isowt2 = 2000*isowt2/bigmax; % Multiply adc by 100 as Andy does adc = 100*adc; for i=1:slices testb(:,:,i)=medfilt2(aniso(:,:,i), [21 21]); end bigmax = max(max(max(testb))); aniso = 2000*aniso/bigmax; aniso(aniso>5000)=0; % Write to disk name=strcat(newname,'.ref'); fid2=fopen(name,'wb'); fwrite(fid2, ref,'uint16'); fclose(fid2); name=strcat(newname,'.aniso'); fid2=fopen(name,'wb'); fwrite(fid2, aniso,'uint16'); fclose(fid2); name=strcat(newname,'.BY'); fid2=fopen(name,'wb'); fwrite(fid2, BY,'uint16'); fclose(fid2); name=strcat(newname,'.BX'); fid2=fopen(name,'wb'); fwrite(fid2, BX,'uint16'); fclose(fid2);