clear all close all im_direction = 'vertical'; % Setup path and file naming convention fileFolder = 'C:\Documents and Settings\Carestream Health\My Documents\Backup'; numDelims = 1; delim = ['_']; % Generate cell array of image names collected at each phase V1_dirOutput = dir(fullfile(fileFolder,'V1_*.*')); V2_dirOutput = dir(fullfile(fileFolder,'V2_*.*')); V3_dirOutput = dir(fullfile(fileFolder,'V3_*.*')); OG_dirOutput = dir(fullfile(fileFolder,'OG_*.*')); % Navigate to directory where images are cd('C:\Documents and Settings\Carestream Health\My Documents\Backup'); % Create variables (cell arrays) for each image filename V1_fileNames = {V1_dirOutput.name}'; V2_fileNames = {V2_dirOutput.name}'; V3_fileNames = {V3_dirOutput.name}'; OG_fileNames = {OG_dirOutput.name}'; % Check the number of files collected at each phase V1_numFiles = size(V1_fileNames,1); V2_numFiles = size(V2_fileNames,1); V3_numFiles = size(V3_fileNames,1); OG_numFiles = size(OG_fileNames,1); % Ensure the the number of files collected at each phase are equal if (isequal(V1_numFiles,V2_numFiles,V3_numFiles,OG_numFiles)==1) disp('Processing Phase 1 images:'); disp([V1_fileNames]); disp('Processing Phase 2 images:'); disp([V2_fileNames]); disp('Processing Phase 3 images:'); disp([V3_fileNames]); disp('Processing OG images:'); disp([OG_fileNames]); m = V1_numFiles; clear('V1_numFiles', 'V2_numFiles', 'V3_numFiles', 'OG_numFiles'); else disp('Error: Number of files processed for each phase and reconstruction must be equal'); end % Create a matrix of the file list index and the delimiters V1_fileNums=[ [1:m]' zeros(m,1) ]; V2_fileNums=[ [1:m]' zeros(m,1) ]; V3_fileNums=[ [1:m]' zeros(m,1) ]; OG_fileNums=[ [1:m]' zeros(m,1) ]; % Cycles through list of files to determine number of iterations collected for i=1:m V1_rem=V1_dirOutput(i).name; V2_rem=V2_dirOutput(i).name; V3_rem=V3_dirOutput(i).name; OG_rem=OG_dirOutput(i).name; [V1_token,V1_rem] = strtok(V1_rem,delim); [V2_token,V2_rem] = strtok(V2_rem,delim); [V3_token,V3_rem] = strtok(V3_rem,delim); [OG_token,OG_rem] = strtok(OG_rem,delim); V1_fileNums(i,3) = str2num(V1_rem(1,2:end-4)); V2_fileNums(i,3) = str2num(V2_rem(1,2:end-4)); V3_fileNums(i,3) = str2num(V3_rem(1,2:end-4)); OG_fileNums(i,3) = str2num(OG_rem(1,2:end-4)); end clear V1_rem V2_rem V3_rem OG_rem V1_token V2_token V3_token OG_token; % % Sort the matrix by rows V1_fileNums = sortrows(V1_fileNums,2); V2_fileNums = sortrows(V2_fileNums,2); V3_fileNums = sortrows(V3_fileNums,2); OG_fileNums = sortrows(OG_fileNums,2); % Show the filenames in the new order V1_dirOutput(V1_fileNums(:,1)).name; V2_dirOutput(V2_fileNums(:,1)).name; V3_dirOutput(V3_fileNums(:,1)).name; OG_dirOutput(OG_fileNums(:,1)).name; % Create a cell array of the filenames in new order for i=1:m V1{i,1} = V1_dirOutput(V1_fileNums(i,1)).name; V2{i,1} = V2_dirOutput(V2_fileNums(i,1)).name; V3{i,1} = V3_dirOutput(V3_fileNums(i,1)).name; OG{i,1} = OG_dirOutput(OG_fileNums(i,1)).name; end clear V1_dirOutput V1_fileNums V2_dirOutput V2_fileNums V3_dirOutput V3_fileNums OG_dirOutput OG_fileNums; % Define name field for raw images fields = ('name'); % Convert files names for Phae 1, 2, and 3 to structures V1 = cell2struct(V1,fields,m); V2 = cell2struct(V2,fields,m); V3 = cell2struct(V3,fields,m); OG = cell2struct(OG,fields,m); for i = 1:m % Raw image read I1_temp = imread(V1(i).name); I2_temp = imread(V2(i).name); I3_temp = imread(V3(i).name); OG_temp = imread(OG(i).name); % Convert to double I1 = double(I1_temp); I2 = double(I2_temp); I3 = double(I3_temp); OG_I = double(OG_temp); clear I1_temp I2_temp I3_temp OG_temp; % Check that image sizes are equal [r1,c1] = size(I1); [r2,c2] = size(I2); [r3,c3] = size(I3); [rog,cog] = size(OG_I); if isequal(r1,r2,r3,rog) && isequal(c1,c2,c3,cog) clear r2 r3 rog c2 c3 cog; else errordlg('Image sizes are not equal','Check Image Sizes'); end % Determine increments for image line sampling if isequal(im_direction,'horizontal') inc = floor(c1/6); start_col = 6; end_col = c1; cols = start_col:inc:end_col; % Average together 6 columns of the images L1_avg = mean(I1(:,cols),2); L2_avg = mean(I2(:,cols),2); L3_avg = mean(I3(:,cols),2); OG_L_avg = mean(OG_I(:,cols),2); % Select 1 col of each image for comparison of single & avg lines L1_single = I1(:,round(c1/2)); L2_single = I2(:,round(c1/2)); L3_single = I3(:,round(c1/2)); OG_L_single = OG_I(:,round(c1/2)); N = r1; else inc = floor(r1/6); start_line = 6; end_line = r1; rows = start_line:inc:end_line; % Average together 6 lines of the images L1_avg = mean(I1(rows,:),1); L2_avg = mean(I2(rows,:),1); L3_avg = mean(I3(rows,:),1); OG_L_avg = mean(OG_I(rows,:),1); % Select 1 line of each image for comparison of single & avg lines L1_single = I1(round(r1/2),:); L2_single = I2(round(r1/2),:); L3_single = I3(round(r1/2),:); OG_L_single = OG_I(round(r1/2),:); N = c1; end clear inc start_line end_line rows I1 I2 I3; clear inc start_col end_col cols; % Force image lines to be column vectors L1_avg = L1_avg(:); L2_avg = L2_avg(:); L3_avg = L3_avg(:); OG_L_avg = OG_L_avg(:); L1_single = L1_single(:); L2_single = L2_single(:); L3_single = L3_single(:); OG_L_single = OG_L_single(:); % Initialize misc. variables for FFT processing Ts=.001; time=[0:N-1]*Ts; werror=2*pi*1e-6; time=time(:)'; Nt=round((max(time)-min(time)+Ts)/Ts); tvecti=round(time/Ts); tvecti=tvecti-min(tvecti)+1; maxCyc=30; % FFT computations Raw 1 I1_avg =zeros(Nt,1); I2_avg =zeros(Nt,1); I3_avg =zeros(Nt,1); OG_I_avg =zeros(Nt,1); I1_single =zeros(Nt,1); I2_single =zeros(Nt,1); I3_single =zeros(Nt,1); OG_I_single =zeros(Nt,1); % Pad the averaged image lines for the FFT I1_avg(tvecti)= L1_avg; I2_avg(tvecti)= L2_avg; I3_avg(tvecti)= L3_avg; OG_I_avg(tvecti)= OG_L_avg; % Pad the single image lines for the FFT I1_single(tvecti)= L1_single; I2_single(tvecti)= L2_single; I3_single(tvecti)= L3_single; OG_I_single(tvecti)= OG_L_single; % Perform 1D fft on averaged image lines I1_avg_fft = fft(I1_avg); I2_avg_fft = fft(I2_avg); I3_avg_fft = fft(I3_avg); OG_avg_fft = fft(OG_I_avg); % Perform 1D fft on single image lines I1_single_fft = fft(I1_single); I2_single_fft = fft(I2_single); I3_single_fft = fft(I3_single); OG_single_fft = fft(OG_I_single); % Compute FFT magnitudes I1_avg_ps = abs(I1_avg_fft); I2_avg_ps = abs(I2_avg_fft); I3_avg_ps = abs(I3_avg_fft); OG_avg_ps = abs(OG_avg_fft); % Compute FFT magnitudes I1_single_ps = abs(I1_single_fft); I2_single_ps = abs(I2_single_fft); I3_single_ps = abs(I3_single_fft); OG_single_ps = abs(OG_single_fft); % Find 1st Harmonic bin index (frequency) [Mfft_avg_1,w_avg_1] = max(I1_avg_ps(2:round(Nt/2))); [Mfft_avg_2,w_avg_2] = max(I2_avg_ps(2:round(Nt/2))); [Mfft_avg_3,w_avg_3] = max(I3_avg_ps(2:round(Nt/2))); [Mfft_avg_4,w_avg_OG] = max(OG_avg_ps(2:round(Nt/2))); % Find 1st Harmonic bin index (frequency) [Mfft_single_1,w_single_1] = max(I1_single_ps(2:round(Nt/2))); [Mfft_single_2,w_single_2] = max(I2_single_ps(2:round(Nt/2))); [Mfft_single_3,w_single_3] = max(I3_single_ps(2:round(Nt/2))); [Mfft_single_4,w_single_OG] = max(OG_single_ps(2:round(Nt/2))); % Translate the frequency into higher resolution sampling space w_avg_1 = 2*pi*w_avg_1/Nt/Ts; w_avg_2 = 2*pi*w_avg_2/Nt/Ts; w_avg_3 = 2*pi*w_avg_3/Nt/Ts; w_avg_4 = 2*pi*w_avg_OG/Nt/Ts; % Translate the frequency into higher resolution sampling space w_single_1 = 2*pi*w_single_1/Nt/Ts; w_single_2 = 2*pi*w_single_2/Nt/Ts; w_single_3 = 2*pi*w_single_3/Nt/Ts; w_single_4 = 2*pi*w_single_OG/Nt/Ts; % Processing for averaged % lines%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize stucture variables to hold info for each acquisition trial for j = 1:4 S(i,j).bin = zeros(3,1); S(i,j).bin(1) = 0; T(i,j).bin = zeros(3,1); T(i,j).bin(1) = 0; end % Set up structure needed to solve for "true" w D0_avg_1 = [cos(w_avg_1*time);sin(w_avg_1*time);ones(1,N)]'; D0_avg_2 = [cos(w_avg_2*time);sin(w_avg_2*time);ones(1,N)]'; D0_avg_3 = [cos(w_avg_3*time);sin(w_avg_3*time);ones(1,N)]'; D0_avg_OG = [cos(w_avg_4*time);sin(w_avg_4*time);ones(1,N)]'; % Set up structure needed to solve for "true" w D0_single_1 = [cos(w_single_1*time);sin(w_single_1*time);ones(1,N)]'; D0_single_2 = [cos(w_single_2*time);sin(w_single_2*time);ones(1,N)]'; D0_single_3 = [cos(w_single_3*time);sin(w_single_3*time);ones(1,N)]'; D0_single_OG = [cos(w_single_4*time);sin(w_single_4*time);ones(1,N)]'; % Solve for sine wave image components first iteration x0_avg_1 = D0_avg_1\L1_avg; x0_avg_2 = D0_avg_2\L2_avg; x0_avg_3 = D0_avg_3\L3_avg; x0_avg_4 = D0_avg_OG\OG_L_avg; % Solve for sine wave image components first iteration x0_single_1 = D0_single_1\L1_single; x0_single_2 = D0_single_2\L2_single; x0_single_3 = D0_single_3\L3_single; x0_single_4 = D0_single_OG\OG_L_single; x0_avg_1(4)=0; x0_avg_2(4)=0; x0_avg_3(4)=0; x0_avg_4(4)=0; x0_single_1(4)=0; x0_single_2(4)=0; x0_single_3(4)=0; x0_single_4(4)=0; icyc=0; % RAW 1: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_avg_1=[cos(w_avg_1*time);sin(w_avg_1*time);ones(1,N);-x0_avg_1(1)*time.*sin(w_avg_1*time)+x0_avg_1(2)*time.*cos(w_avg_1*time)]'; x0_avg_1=D0_avg_1\L1_avg; w_avg_1=w_avg_1+x0_avg_1(4); if abs(x0_avg_1(4))=maxCyc break; end end % RAW 1: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_single_1=[cos(w_single_1*time);sin(w_single_1*time);ones(1,N);-x0_single_1(1)*time.*sin(w_single_1*time)+x0_single_1(2)*time.*cos(w_single_1*time)]'; x0_single_1=D0_single_1\L1_single; w_single_1=w_single_1+x0_single_1(4); if abs(x0_single_1(4))=maxCyc break; end end % RAW 2: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_avg_2=[cos(w_avg_2*time);sin(w_avg_2*time);ones(1,N);-x0_avg_2(1)*time.*sin(w_avg_2*time)+x0_avg_2(2)*time.*cos(w_avg_2*time)]'; x0_avg_2=D0_avg_2\L2_avg; w_avg_2=w_avg_2+x0_avg_2(4); if abs(x0_avg_2(4))=maxCyc break; end end % RAW 2: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_single_2=[cos(w_single_2*time);sin(w_single_2*time);ones(1,N);-x0_single_2(1)*time.*sin(w_single_2*time)+x0_single_2(2)*time.*cos(w_single_2*time)]'; x0_single_2=D0_single_2\L2_single; w_single_2=w_single_2+x0_single_2(4); if abs(x0_single_2(4))=maxCyc break; end end % RAW 3: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_avg_3=[cos(w_avg_3*time);sin(w_avg_3*time);ones(1,N);-x0_avg_3(1)*time.*sin(w_avg_3*time)+x0_avg_3(2)*time.*cos(w_avg_3*time)]'; x0_avg_3=D0_avg_3\L3_avg; w_avg_3=w_avg_3+x0_avg_3(4); if abs(x0_avg_3(4))=maxCyc break; end end % RAW 3: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_single_3=[cos(w_single_3*time);sin(w_single_3*time);ones(1,N);-x0_single_3(1)*time.*sin(w_single_3*time)+x0_single_3(2)*time.*cos(w_single_3*time)]'; x0_single_3=D0_single_3\L3_single; w_single_3=w_single_3+x0_single_3(4); if abs(x0_single_3(4))=maxCyc break; end end % OG: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_avg_4=[cos(w_avg_4*time);sin(w_avg_4*time);ones(1,N);-x0_avg_4(1)*time.*sin(w_avg_4*time)+x0_avg_4(2)*time.*cos(w_avg_4*time)]'; x0_avg_4=D0_avg_4\OG_L_avg; w_avg_4=w_avg_4+x0_avg_4(4); if abs(x0_avg_4(4))=maxCyc break; end end % OG: Solve for sine wave image components second iteration while 1 icyc=icyc+1; D0_single_4=[cos(w_single_4*time);sin(w_single_4*time);ones(1,N);-x0_single_4(1)*time.*sin(w_single_4*time)+x0_single_4(2)*time.*cos(w_single_4*time)]'; x0_single_4=D0_single_4\OG_L_single; w_single_4=w_single_4+x0_single_4(4); if abs(x0_single_4(4))=maxCyc break; end end % RAW 1: Unwrapped Phase angle in radians of fitted sine wave if x0_avg_1(1)>0 S(i,1).phi_avg = atan(-x0_avg_1(2)/x0_avg_1(1)); elseif x0_avg_1(1)<0 S(i,1).phi_avg = atan(-x0_avg_1(2)/x0_avg_1(1))+ pi; else if x0_avg_1(2)<0 S(i,1).phi_avg=pi/2; else S(i,1).phi_avg=3*pi/2; end end % Convert to degrees S(i,1).phi_avg=S(i,1).phi_avg*(180/pi); % RAW 1: Unwrapped Phase angle in radians of fitted sine wave if x0_single_1(1)>0 S(i,1).phi_single = atan(-x0_single_1(2)/x0_single_1(1)); elseif x0_single_1(1)<0 S(i,1).phi_single = atan(-x0_single_1(2)/x0_single_1(1))+ pi; else if x0_single_1(2)<0 S(i,1).phi_single=pi/2; else S(i,1).phi_single=3*pi/2; end end % Convert to degrees S(i,1).phi_single=S(i,1).phi_single*(180/pi); % RAW 2: Unwrapped Phase angle in radians of fitted sine wave if x0_avg_2(1)>0 S(i,2).phi_avg = atan(-x0_avg_2(2)/x0_avg_2(1)); elseif x0_avg_2(1)<0 S(i,2).phi_avg = atan(-x0_avg_2(2)/x0_avg_2(1))+ pi; else if x0_avg_2(2)<0 S(i,2).phi_avg=pi/2; else S(i,2).phi_avg=3*pi/2; end end % Convert to degrees S(i,2).phi_avg=S(i,2).phi_avg*(180/pi); % RAW 2: Unwrapped Phase angle in radians of fitted sine wave if x0_single_2(1)>0 S(i,2).phi_single = atan(-x0_single_2(2)/x0_single_2(1)); elseif x0_single_2(1)<0 S(i,2).phi_single = atan(-x0_single_2(2)/x0_single_2(1))+ pi; else if x0_single_2(2)<0 S(i,2).phi_single=pi/2; else S(i,2).phi_single=3*pi/2; end end % Convert to degrees S(i,2).phi_single=S(i,2).phi_single*(180/pi); % RAW 3: Unwrapped Phase angle in radians of fitted sine wave if x0_avg_3(1)>0 S(i,3).phi_avg = atan(-x0_avg_3(2)/x0_avg_3(1)); elseif x0_avg_3(1)<0 S(i,3).phi_avg = atan(-x0_avg_3(2)/x0_avg_3(1))+ pi; else if x0_avg_3(2)<0 S(i,3).phi_avg=pi/2; else S(i,3).phi_avg=3*pi/2; end end % Convert to degrees S(i,3).phi_avg=S(i,3).phi_avg*(180/pi); % RAW 3: Unwrapped Phase angle in radians of fitted sine wave if x0_single_3(1)>0 S(i,3).phi_single = atan(-x0_single_3(2)/x0_single_3(1)); elseif x0_single_3(1)<0 S(i,3).phi_single = atan(-x0_single_3(2)/x0_single_3(1))+ pi; else if x0_single_3(2)<0 S(i,3).phi_single=pi/2; else S(i,3).phi_single=3*pi/2; end end % Convert to degrees S(i,3).phi_single=S(i,3).phi_single*(180/pi); % OG: Unwrapped Phase angle in radians of fitted sine wave if x0_avg_4(1)>0 S(i,4).phi_avg = atan(-x0_avg_4(2)/x0_avg_4(1)); elseif x0_avg_4(1)<0 S(i,4).phi_avg = atan(-x0_avg_4(2)/x0_avg_4(1))+ pi; else if x0_avg_4(2)<0 S(i,4).phi_avg=pi/2; else S(i,4).phi_avg=3*pi/2; end end % Convert to degrees S(i,4).phi_avg=S(i,4).phi_avg*(180/pi); % OG: Unwrapped Phase angle in radians of fitted sine wave if x0_single_4(1)>0 S(i,4).phi_single = atan(-x0_single_4(2)/x0_single_4(1)); elseif x0_single_4(1)<0 S(i,4).phi_single = atan(-x0_single_4(2)/x0_single_4(1))+ pi; else if x0_single_4(2)<0 S(i,4).phi_single=pi/2; else S(i,4).phi_single=3*pi/2; end end % Convert to degrees S(i,4).phi_single=S(i,4).phi_single*(180/pi); % Real component of fitted sine wave S(i,1).real_avg = x0_avg_1(1); S(i,2).real_avg = x0_avg_2(1); S(i,3).real_avg = x0_avg_3(1); S(i,4).real_avg = x0_avg_4(1); % Imag component of fitted sine wave S(i,1).imag_avg = x0_avg_1(2); S(i,2).imag_avg = x0_avg_2(2); S(i,3).imag_avg = x0_avg_3(2); S(i,4).imag_avg = x0_avg_4(2); % DC component of fitted sine wave (1/pixel value) S(i,1).DC_avg=x0_avg_1(3); S(i,2).DC_avg=x0_avg_2(3); S(i,3).DC_avg=x0_avg_3(3); S(i,4).DC_avg=x0_avg_4(3); % Amplitude of fitted sine wave S(i,1).A_avg=(sqrt(x0_avg_1(1)^2+x0_avg_1(2)^2)); S(i,2).A_avg=(sqrt(x0_avg_2(1)^2+x0_avg_2(2)^2)); S(i,3).A_avg=(sqrt(x0_avg_3(1)^2+x0_avg_3(2)^2)); S(i,4).A_avg=(sqrt(x0_avg_4(1)^2+x0_avg_4(2)^2)); % Frequency of sine wave in Hz S(i,1).f_avg = w_avg_1/2/pi; S(i,2).f_avg = w_avg_2/2/pi; S(i,3).f_avg = w_avg_3/2/pi; S(i,4).f_avg = w_avg_4/2/pi; % FFT of image line shifted so DC is at center S(i,1).fft_avg = fftshift(I1_avg_ps*Ts); S(i,2).fft_avg = fftshift(I2_avg_ps*Ts); S(i,3).fft_avg = fftshift(I3_avg_ps*Ts); S(i,4).fft_avg = fftshift(OG_avg_ps*Ts); figure(1) plot(S(i,1).fft_avg,'-b*'); hold on; plot(S(i,2).fft_avg,'-m*'); plot(S(i,3).fft_avg,'-g*'); plot(S(i,4).fft_avg,'-ks'); hold off; % Real component of fitted sine wave S(i,1).real_single = x0_single_1(1); S(i,2).real_single = x0_single_2(1); S(i,3).real_single = x0_single_3(1); S(i,4).real_single = x0_single_4(1); % Imag component of fitted sine wave S(i,1).imag_single = x0_single_1(2); S(i,2).imag_single = x0_single_2(2); S(i,3).imag_single = x0_single_3(2); S(i,4).imag_single = x0_single_4(2); % DC component of fitted sine wave (1/pixel value) S(i,1).DC_single=x0_single_1(3); S(i,2).DC_single=x0_single_2(3); S(i,3).DC_single=x0_single_3(3); S(i,4).DC_single=x0_single_4(3); % Amplitude of fitted sine wave S(i,1).A_single=(sqrt(x0_single_1(1)^2+x0_single_1(2)^2)); S(i,2).A_single=(sqrt(x0_single_2(1)^2+x0_single_2(2)^2)); S(i,3).A_single=(sqrt(x0_single_3(1)^2+x0_single_3(2)^2)); S(i,4).A_single=(sqrt(x0_single_4(1)^2+x0_single_4(2)^2)); % Frequency of sine wave in Hz S(i,1).f_single = w_single_1/2/pi; S(i,2).f_single = w_single_2/2/pi; S(i,3).f_single = w_single_3/2/pi; S(i,4).f_single = w_single_4/2/pi; % FFT of image line shifted so DC is at center S(i,1).fft_single = fftshift(I1_single_ps*Ts); S(i,2).fft_single = fftshift(I2_single_ps*Ts); S(i,3).fft_single = fftshift(I3_single_ps*Ts); S(i,4).fft_single = fftshift(OG_single_ps*Ts); figure(2) plot(S(i,1).fft_single,'-b*'); hold on; plot(S(i,2).fft_single,'-m*'); plot(S(i,3).fft_single,'-g*'); plot(S(i,4).fft_single,'-ks'); hold off; for k = 0:3 % Raw 1 & OG bin_lo(k+1) = round(N/2 + k*S(i,1).f_avg)-5; bin_hi(k+1) = round(N/2 + k*S(i,1).f_avg)+5; range = (bin_lo(k+1)+1):(bin_hi(k+1)+1); [S(i,1).h_avg(k+1),bin_temp] = max(S(i,1).fft_avg(range,1)); S(i,1).bin_avg(k+1) = range(bin_temp); [S(i,4).h_avg(k+1),bin_temp] = max(S(i,4).fft_avg(range,1)); S(i,4).bin_avg(k+1)= range(bin_temp); clear 'bin_hi' 'bin_lo' 'range' 'bin_temp'; % Raw 2 bin_lo(k+1) = round(N/2 + k*S(i,2).f_avg)-5; bin_hi(k+1) = round(N/2 + k*S(i,2).f_avg)+5; range = (bin_lo(k+1)+1):(bin_hi(k+1)+1); [S(i,2).h_avg(k+1),bin_temp] = max(S(i,2).fft_avg(range,1)); S(i,2).bin_avg(k+1) = range(bin_temp); clear 'bin_hi' 'bin_lo' 'range' 'bin_temp'; % Raw 3 bin_lo(k+1) = round(N/2 + k*S(i,3).f_avg)-5; bin_hi(k+1) = round(N/2 + k*S(i,3).f_avg)+5; range = (bin_lo(k+1)+1):(bin_hi(k+1)+1); [S(i,3).h_avg(k+1),bin_temp] = max(S(i,3).fft_avg(range,1)); S(i,3).bin_avg(k+1) = range(bin_temp); clear 'bin_hi' 'bin_lo' 'range' 'bin_temp'; end for j = 1:4 S(i,j).merit_avg = 100.*(S(i,4).h_avg(3)/S(i,1).h_avg(2)); end for k = 0:3 % Raw 1 & OG bin_lo(k+1) = round(N/2 + k*S(i,1).f_single)-5; bin_hi(k+1) = round(N/2 + k*S(i,1).f_single)+5; range = (bin_lo(k+1)+1):(bin_hi(k+1)+1); [S(i,1).h_single(k+1),bin_temp] = max(S(i,1).fft_single(range,1)); S(i,1).bin_single(k+1) = range(bin_temp); [S(i,4).h_single(k+1),bin_temp] = max(S(i,4).fft_single(range,1)); S(i,4).bin_single(k+1)= range(bin_temp); clear 'bin_hi' 'bin_lo' 'range' 'bin_temp'; % Raw 2 bin_lo(k+1) = round(N/2 + k*S(i,2).f_single)-5; bin_hi(k+1) = round(N/2 + k*S(i,2).f_single)+5; range = (bin_lo(k+1)+1):(bin_hi(k+1)+1); [S(i,2).h_single(k+1),bin_temp] = max(S(i,2).fft_single(range,1)); S(i,2).bin_single(k+1) = range(bin_temp); clear 'bin_hi' 'bin_lo' 'range' 'bin_temp'; % Raw 3 bin_lo(k+1) = round(N/2 + k*S(i,3).f_single)-5; bin_hi(k+1) = round(N/2 + k*S(i,3).f_single)+5; range = (bin_lo(k+1)+1):(bin_hi(k+1)+1); [S(i,3).h_single(k+1),bin_temp] = max(S(i,3).fft_single(range,1)); S(i,3).bin_single(k+1) = range(bin_temp); clear 'bin_hi' 'bin_lo' 'range' 'bin_temp'; end for j = 1:4 S(i,j).merit_single = 100.*(S(i,4).h_single(3)/S(i,1).h_single(2)); end end BZ = [4.751 4.759 4.768 4.772 4.788 4.787... 4.749 4.755 4.770 4.766 4.785 4.795... 4.752 4.757 4.779 4.771 4.789 4.791... 4.749 4.683 4.776 4.772 4.790 4.782... 4.756 4.767 4.772 4.769 4.788 4.780... 4.756 4.764 4.770 4.772 4.784 4.791]; MM = 100.*[0.1506 0.1485 0.1467 0.1478 0.1470 0.1450... 0.1483 0.1486 0.1507 0.1493 0.1488 0.1480... 0.1451 0.1516 0.1451 0.1479 0.1479 0.1473... 0.1479 0.1483 0.1483 0.1483 0.1492 0.1452... 0.1454 0.1468 0.1498 0.1501 0.1459 0.1474... 0.1452 0.1479 0.1466 0.1516 0.1473 0.1469]; for q = 1:4 %BJM Changed from 36 to 45 for j = 1:36 S(j,q).meritBZ = BZ(1,j); S(j,q).meritMM = MM(1,j); end end for i = 1:m P1_avg(i) = S(i,1).phi_avg; P2_avg(i) = S(i,2).phi_avg; P3_avg(i) = S(i,3).phi_avg; merit_P_avg(i) = S(i,1).merit_avg; % merit_PBZ(i) = S(i,1).meritBZ; % merit_PMM(i) = S(i,1).meritMM; P1_single(i) = S(i,1).phi_single; P2_single(i) = S(i,2).phi_single; P3_single(i) = S(i,3).phi_single; merit_P_single(i) = S(i,1).merit_single; end P2_min_avg = min(P2_avg); P2_max_avg = max(P2_avg); P3_min_avg = min(P3_avg); P3_max_avg = max(P3_avg); P2_min_single = min(P2_single); P2_max_single = max(P2_single); P3_min_single = min(P3_single); P3_max_single = max(P3_single); phase_matrix_avg = [P1_avg;P2_avg;P3_avg;merit_P_avg]; phase_matrix_single = [P1_single;P2_single;P3_single;merit_P_single]; saveFolder = 'C:\Documents and Settings\Carestream Health\My Documents\Backup'; savefile = fullfile(saveFolder,'trial_1.mat'); %BJM commented out below line to eliminate BZ and MM %save(savefile,'P1_avg','P2_avg','P3_avg','merit_P_avg','P1_single','P2_single','P3_single','merit_P_single','merit_PBZ','merit_PMM'); save(savefile,'P1_avg','P2_avg','P3_avg','merit_P_avg','P1_single','P2_single','P3_single','merit_P_single'); V2_coarse = 44.87:.8:48.07; V3_coarse = [70.58 70.78 71.58 72.38 72.58] [V2,V3] = meshgrid(V2_coarse,V3_coarse); [r_V2,c_V2] = size(V2_coarse); [r_V3,c_V3] = size(V3_coarse); k = 1; for i = 1:c_V3 for j = 1:c_V2 merit_score_avg(i,j) = S(k,1).merit_avg; merit_score_single(i,j) = S(k,1).merit_single; merit_score_BZ(i,j) = S(k,1).meritBZ; merit_score_MM(i,j) = S(k,1).meritMM; k = k+1; end end figure(3) surf(V3,V2,merit_score_avg); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); title('Merit Map SAM Technique Avg 6 Lines'); legend('6 Line Avg'); colormap('hot'); hold off; figure(4) surf(V3,V2,merit_score_single); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); title('Merit Map SAM Technique Single Line'); legend('Single Line'); colormap('hot'); hold off; figure(5) surf(V3,V2,merit_score_avg); hold on; surf(V3,V2,merit_score_single); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); title('Merit Map SAM Technique Avg Compared to Single Line'); legend('6 Line Avg','Single Line'); colormap('hot'); hold off; figure(6) surf(V3,V2,merit_score_BZ); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); legend('BZ Excel'); title('Merit Map SAM Technique, BZ Excel Data'); colormap('hot'); hold off figure(7) surf(V3,V2,merit_score_avg); hold on; surf(V3,V2,merit_score_single); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); title('Merit Map SAM Technique, Avg Compared to Single Line and BZ Excel Data'); hold on; surf(V3,V2,merit_score_BZ); legend('6 Line Avg','Single Line','BZ Excel'); colormap('hot'); hold off figure(8) surf(V3,V2,merit_score_MM); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); legend('MetaMorph'); title('Merit Map SAM Technique, MetaMorph Merit Data'); colormap('hot'); hold off figure(9) surf(V3,V2,merit_score_avg); hold on; surf(V3,V2,merit_score_single); surf(V3,V2,merit_score_BZ); surf(V3,V2,merit_score_MM); xlabel('V3 Range (volts)'); ylabel('V2 Range (volts)'); zlabel('Merit Score (%)'); legend('6 Line Avg','Single Line','BZ Excel','MetaMorph'); title('Merit Map SAM Technique, Merit Data Comparison'); colormap('hot'); hold off % % xlin = linspace(P3_min,P3_max,.5); % ylin = linspace(P2_min,P2_max,.5); % [X,Y] = meshgrid(ylin,xlin); % Z = griddata(V3,V2,merit_score,X,Y,'cubic'); % figure(5) % mesh(X,Y,Z); % hold on; % axis tight; % plot3(P3,P2,merit_P,'.','MarkerSize',15); % hold off % % %