diff --git a/GLMdenoisedata.m b/GLMdenoisedata.m
index 1b2f1cf..e4324fa 100644
--- a/GLMdenoisedata.m
+++ b/GLMdenoisedata.m
@@ -172,6 +172,10 @@
 %     column vectors (XYZ x 1) into a more palatable form. Default is to
 %     do nothing special: @(vals) vals. Note that we do not save this input
 %     to results.inputs.opt in order to avoid weird .mat file-saving problems.
+%   <reconmask> (optional) is a X x Y x Z mask that will be used to turn the
+%     (presumably) 3D data into a data x time vector (save RAM, time, sanity)
+%     also used to reconstruct the data back into its 3D state for figures.
+%   <wantpcmap> (optional) is whether to compute and save PC maps. Default: 0.
 % <figuredir> (optional) is a directory to which to write figures.  (If the
 %   directory does not exist, we create it; if the directory already exists,
 %   we delete its contents so we can start afresh.)  If [], no figures are
@@ -364,7 +368,7 @@
 %   that the y-axis is now converted into the equivalent amount of data gained.  
 %   For example, if SNR is boosted from 4 to 8, this is equivalent to having 
 %   obtained 300% more data than was actually obtained.
-% - "PCmap/PCmap_runN_numO.png" shows the estimated weights for the Oth PC
+% - "PCmap/PCmap_runN_numO.png" (if wantpcmap==1) shows the estimated weights for the Oth PC
 %   for run N.  The range is -A to A where A is the 99th percentile of the
 %   absolute value of all weights across all runs.  The colormap proceeds
 %   from blue (negative) to black (0) to red (positive).
@@ -382,6 +386,13 @@
 % times at which data are actually sampled.
 %
 % History:
+% - 2025/08/20: Having an issue with PCMap plotting, off and on. Adding option 
+%               to disable until I can sort that out. 
+% - 2018/07/20: Added plots of PC timecourse below pcweight images.  
+%               Added a new struct for opt, 'reconmask'. This will be  
+%               used to vectorize the data. Suggestion is whole brain mask
+%               (including whitematter/CSF) as a binary map with the exact same
+%               dimensions as data. 3D images will be reconstructed for figure creation
 % - 2016/09/02: to avoid weird .mat file saving issues, do not save inputs.opt.drawfunction
 % - 2016/04/15: add opt.drawfunction
 % - 2014/08/01: add opt.wantparametric input (which enables parametric GLM fits).
@@ -470,6 +481,39 @@
   fprintf('*** GLMdenoisedata: WARNING: we checked the first run and found some non-finite values (e.g. NaN, Inf). unexpected results may occur due to non-finite values. please fix and re-run GLMdenoisedata. ***\n');
 end
 
+% if opt.reconmask is provided, vectorize and remove zeros
+if isfield(opt,'reconmask') 
+    numvoxels = size(reshape(opt.reconmask, [],1),1);
+    fprintf('*** GLMdenoisedata: vectorizing data.');
+    mask_ind = opt.reconmask(:,:,:,1) ~= 0; % generate the indices where data live
+    for i = 1:size(data,2)
+        t_dim = size(data{1,i},4);
+        data{1,i} = reshape(data{1,i}, [], t_dim);
+        data{1,i} = data{1,i}(mask_ind, :); %generates a vector containing values within mask.
+        new_vox_count = size(data{1,i},1);
+        fprintf('.')
+    end
+
+    if isfield(opt, 'hrffitmask') %convert the provided data masking volumes to equivilent vectors
+    	opt.hrffitmask = reshape(opt.hrffitmask, [], 1)
+    	opt.hrffitmask = opt.hrffitmask(mask_ind, :);
+    end
+
+    if isfield(opt, 'brainexclude')
+    	opt.brainexclude = reshape(opt.brainexclude, [] ,1)
+    	opt.brainexclude = opt.brainexclude(mask_ind, :);
+    end
+
+    if isfield(opt, 'pcR2cutoffmask')
+    	opt.pcR2cutoffmask = reshape(opt.pcR2cutoffmask, [] ,1)
+    	opt.pcR2cutoffmask = opt.pcR2cutoffmask(mask_ind, :);
+    end
+
+    vec_reduced = (new_vox_count/numvoxels) *100;
+    fprintf('\nData reduced to %.2f%% of original size using provided mask \n', vec_reduced);
+end
+
+
 % calc
 numruns = length(design);
 dataclass = class(data{1});  % will always be 'single'
@@ -555,6 +599,9 @@
 if ~isfield(opt,'wantsanityfigures') || isempty(opt.wantsanityfigures)
   opt.wantsanityfigures = 1;
 end
+if ~isfield(opt,'wantpcmap') || isempty(opt.wantpcmap)
+  opt.wantpcmap = 0;
+end
 if ~isfield(opt,'drawfunction') || isempty(opt.drawfunction)
   opt.drawfunction = @(vals) vals;
 end
@@ -574,7 +621,9 @@
     assert(rmdir(figuredir,'s'));
   end
   assert(mkdir(figuredir));
-  assert(mkdir(fullfile(figuredir,'PCmap')));
+  if opt.wantpcmap
+    assert(mkdir(fullfile(figuredir,'PCmap')));
+  end
   if opt.wantsanityfigures
     assert(mkdir(fullfile(figuredir,'CheckData')));
   end
@@ -1073,8 +1122,24 @@
     figurewrite('HRF',[],[],figuredir);
   end
   
+   %create the correctly sized mask for reconstructing 3d images from vecs
+    if isfield(opt, 'reconmask')
+        mask_vec = reshape(opt.reconmask,[],1);
+        x_dim = size(opt.reconmask,1);
+        y_dim = size(opt.reconmask,2);
+        z_dim = size(opt.reconmask,3);
+    end
+    
+  
   % write out image showing HRF fit voxels
   if isequal(hrfmodel,'optimize') && ~isempty(results.hrffitvoxels)
+        if isfield(opt, 'reconmask')
+            mask_vec = reshape(opt.reconmask,[],1);
+            hrffit_recon = mask_vec;
+            hrffit_recon(mask_vec ~=0) = results.hrffitvoxels;
+            hrffit_recon = single(reshape(hrffit_recon,x_dim,y_dim,z_dim));
+            results.hrffitvoxels = hrffit_recon;
+        end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(results.hrffitvoxels),[0 1])),gray(256),fullfile(figuredir,'HRFfitvoxels.png'));
   end
 
@@ -1109,20 +1174,44 @@
   end
 
   % write out image showing mean volume (of first run)
+  if isfield(opt, 'reconmask')
+        mean_vol_recon = mask_vec;
+        mean_vol_recon(mask_vec ~=0) = results.meanvol;
+        mean_vol_recon = reshape(mean_vol_recon, x_dim,y_dim,z_dim);
+        results.meanvol=mean_vol_recon;
+  end
   imwrite(uint8(255*makeimagestack(opt.drawfunction(results.meanvol),1)),gray(256),fullfile(figuredir,'MeanVolume.png'));
 
   % write out image showing noise pool
   if ~isempty(results.noisepool)  % in certain degenerate cases, this might happen
+        if isfield(opt, 'reconmask')
+            noisepool_recon = mask_vec;
+            noisepool_recon(mask_vec ~=0) = results.noisepool;
+            noisepool_recon = reshape(noisepool_recon, x_dim,y_dim,z_dim);
+            results.noisepool=noisepool_recon;
+        end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(results.noisepool),[0 1])),gray(256),fullfile(figuredir,'NoisePool.png'));
   end
   
   % write out image showing voxels excluded from noise pool
   if ~isequal(opt.brainexclude,0)
+        if isfield(opt, 'reconmask')
+            temp_recon = mask_vec;
+            temp_recon(mask_vec ~=0) = opt.brainexclude;
+            temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+            opt.brainexclude=temp_recon;
+        end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(opt.brainexclude),[0 1])),gray(256),fullfile(figuredir,'NoiseExclude.png'));
   end
 
   % write out image showing voxel mask for HRF fitting
   if isequal(hrfmodel,'optimize') && ~isequal(opt.hrffitmask,1)
+       if isfield(opt, 'reconmask')
+            temp_recon = mask_vec;
+            temp_recon(mask_vec ~=0) = opt.hrffitmask;
+            temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+            opt.hrffitmask=temp_recon;
+        end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(opt.hrffitmask),[0 1])),gray(256),fullfile(figuredir,'HRFfitmask.png'));
   end
 
@@ -1134,10 +1223,22 @@
 
     % write out image showing voxel mask for PC selection
     if ~isequal(opt.pcR2cutoffmask,1)
+            if isfield(opt, 'reconmask')
+                temp_recon = mask_vec;
+                temp_recon(mask_vec ~=0) = opt.pcR2cutoffmask;
+                temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+                opt.pcR2cutoffmask=temp_recon;
+            end
       imwrite(uint8(255*makeimagestack(opt.drawfunction(opt.pcR2cutoffmask),[0 1])),gray(256),fullfile(figuredir,'PCmask.png'));
     end
   
     % write out image showing the actual voxels used for PC selection
+    if isfield(opt, 'reconmask')
+            temp_recon = mask_vec;
+            temp_recon(mask_vec ~=0) = results.pcvoxels;
+            temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+            results.pcvoxels=temp_recon;
+    end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(results.pcvoxels),[0 1])),gray(256),fullfile(figuredir,'PCvoxels.png'));
 
     % figure out bounds for the R^2 values
@@ -1153,6 +1254,12 @@
     % write out cross-validated R^2 for the various numbers of PCs
     for p=1:size(results.pcR2,dimdata+1)
       temp = subscript(results.pcR2,[repmat({':'},[1 dimdata]) {p}]);
+      if isfield(opt, 'reconmask')
+                temp_recon = mask_vec;
+                temp_recon(mask_vec ~=0) = temp;
+                temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+                temp=temp_recon;
+            end
       feval(imfun,temp,fullfile(figuredir,sprintf('PCcrossvalidation%02d.png',p-1)));
       feval(imfunB,temp,fullfile(figuredir,sprintf('PCcrossvalidationscaled%02d.png',p-1)));
     end
@@ -1160,18 +1267,48 @@
   end
 
   % write out overall R^2 for final model
+  if isfield(opt, 'reconmask')
+        temp_recon = mask_vec;
+        temp_recon(mask_vec ~=0) = results.R2;
+        temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+        results.R2=temp_recon;
+  end
   feval(imfun,results.R2,fullfile(figuredir,sprintf('FinalModel.png')));
 
   % write out R^2 separated by runs for final model
   for p=1:size(results.R2run,dimdata+1)
     temp = subscript(results.R2run,[repmat({':'},[1 dimdata]) {p}]);
+    if isfield(opt, 'reconmask')
+            temp_recon = mask_vec;
+            temp_recon(mask_vec ~=0) = temp;
+            temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+            temp=temp_recon;
+        end
     feval(imfun,temp,fullfile(figuredir,sprintf('FinalModel_run%02d.png',p)));
   end
   
   % write out signal, noise, and SNR
+  if isfield(opt, 'reconmask')
+        temp_recon = mask_vec;
+        temp_recon(mask_vec ~=0) = results.signal;
+        temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+        results.signal=temp_recon;
+  end
   imwrite(uint8(255*makeimagestack(opt.drawfunction(results.signal),[0 prctile(results.signal(:),99)])),hot(256),fullfile(figuredir,'SNRsignal.png'));
   if opt.numboots ~= 0
+                if isfield(opt, 'reconmask')
+                    temp_recon = mask_vec;
+                    temp_recon(mask_vec ~=0) = results.noise;
+                    temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+                    results.noise=temp_recon;
+                end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(results.noise),[0 max(eps,prctile(results.noise(:),99))])),hot(256),fullfile(figuredir,'SNRnoise.png'));
+                if isfield(opt, 'reconmask')
+                    temp_recon = mask_vec;
+                    temp_recon(mask_vec ~=0) = results.SNR;
+                    temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+                    results.SNR=temp_recon;
+                end
     imwrite(uint8(255*makeimagestack(opt.drawfunction(results.SNR),[0 10])),hot(256),fullfile(figuredir,'SNR.png'));
   end
   
@@ -1207,13 +1344,33 @@
   end
   
   % write out maps of pc weights
-  thresh = prctile(abs(results.pcweights(:)),99);
-  for p=1:size(results.pcweights,dimdata+1)
-    for q=1:size(results.pcweights,dimdata+2)
-      temp = subscript(results.pcweights,[repmat({':'},[1 dimdata]) {p} {q}]);
-      imwrite(uint8(255*makeimagestack(opt.drawfunction(temp),[-thresh thresh])),cmapsign(256), ...
-              fullfile(figuredir,'PCmap',sprintf('PCmap_run%02d_num%02d.png',q,p)));
+  if opt.wantpcmap ~= 0
+    thresh = prctile(abs(results.pcweights(:)),99);
+    for p=1:size(results.pcweights,dimdata+1)
+      for q=1:size(results.pcweights,dimdata+2)
+        temp = subscript(results.pcweights,[repmat({':'},[1 dimdata]) {p} {q}]);
+                      if isfield(opt, 'reconmask')
+                          temp_recon = mask_vec;
+                          temp_recon(mask_vec ~=0) = temp;
+                          temp_recon = reshape(temp_recon, x_dim,y_dim,z_dim);
+                          temp=temp_recon;
+                      end
+                      
+                      %Go ahead and create the map as it was originaly done
+                      temp = (255*makeimagestack(opt.drawfunction(temp),[-thresh thresh]));
+                      
+                      %Make sure the x axis is the right length
+                      time_plot_dim = size(results.pcregressors{1,q}(:,p),1); 
+                    
+                      %Create a 3x3 plot space, use the majority for the
+                      %weightmap, only a the last row for weight plot.
+                      subplot(4,3,[1, 2, 3, 4, 5, 6, 7, 8, 9]); imshow(temp,cmapsign(256));
+                      subplot(4,3,[10, 11, 12]); plot(results.pcregressors{1,q}(:,p)); xlim([0, time_plot_dim]);
+                      
+                      %save everything to the right place
+                      label = strcat('PCmap_run',num2str(q), '_num' , num2str(p));
+                      print(fullfile(figuredir, 'PCmap', label), '-dpng');
+      end
     end
-  end
-
+  end % if we want PCmaps or not. 
 end
diff --git a/utilities/imageoverlayer.m b/utilities/imageoverlayer.m
new file mode 100644
index 0000000..b43200d
--- /dev/null
+++ b/utilities/imageoverlayer.m
@@ -0,0 +1,37 @@
+function imageoverlayer(underlay, overlay, thresh)
+    
+    % TODO - make some guess as to if volumes need to be skipped. 
+    skip = 1;
+
+    if size(underlay,3) >1
+
+        ax1 = axes;
+        imagesc(makeimagestack(underlay(:,:,1:skip:end)));
+        colormap(ax1, 'gray')
+        axis('image');
+        hold on
+        % lets mask some things, to see stuff. 
+        alphahelper = overlay(:,:,1:skip:end)>thresh;
+        alphahelper = alphahelper + (overlay(:,:,1:skip:end)<(-1*thresh));
+        ax2 = axes;
+        imagesc(ax2, makeimagestack((overlay(:,:,1:skip:end))), 'AlphaData', makeimagestack(alphahelper))
+        axis('image')
+        colormap(ax2, 'jet')
+        ax2.Visible = 'off';
+        
+    else
+        % its just a slice, get crazy
+        ax1 = axes;
+        imagesc(underlay);
+        colormap(ax1, 'gray')
+        axis('image')
+        hold on
+        % lets mask some things, to see stuff. 
+        alphahelper = overlay>thresh;
+        alphahelper = alphahelper + (overlay<(-1*thresh));
+        ax2 = axes;
+        imagesc(ax2, overlay, 'AlphaData', alphahelper)
+        axis('image')
+        colormap(ax2, 'jet')
+        ax2.Visible = 'off';
+    end