addons commit

Addons/ContrastRatio.m

+function medianCR = ContrastRatio(MatOut)
+
+
+
+%Size of patchs
+SZ = 64;
+SX = 64;
+
+PWTD = CropZeros(MatOut);
+[Nz, Nx] = size(PWTD);
+
+
+PWTD = PWTD(1:floor(Nz/SZ)*SZ,1+Nx-floor(Nx/SX)*SX:end); 
+[Nz, Nx] = size(PWTD);
+
+N11 = SZ*ones(1,Nz/SZ);
+N22 = SX*ones(1,Nx/SX);
+PWTD2 = mat2cell(PWTD,N11,N22);
+PWTD2 = reshape(PWTD2,Nz/SZ*Nx/SX,[]);
+
+%Define reference Patch
+%R21 = PWTD(zz:zz+SZ-1,xx:xx+SX-1);  % fixed position
+patchs=zeros(1,length(PWTD2));
+for i=1:length(PWTD2)
+R11 = PWTD2{i};
+if IsPatchValid(R11)
+patchs(i) = mean(R11(:).^2); 
+else 
+    patchs(i)=0;
+end
+end
+
+patchs(patchs == 0 ) = NaN;
+[A,I] = min(patchs);
+R21 = PWTD2{I};
+RefPower = mean(R21(:).^2);
+
+% Calc CR
+CR=zeros(1,length(PWTD2));
+for i=1:length(PWTD2)
+R11 = PWTD2{i};
+CR(i) = 10*log10(mean(R11(:).^2) / RefPower);  %  = mean(R21(:).^2)
+end
+CR = CR(isfinite(CR));
+CR = CR(CR>0);
+
+medianCR = median(CR,'all');
+
+% show boxplot
+% listis= boxplot(CR, 'Whisker', 2);
+% CRis=get(listis(6)).YData;
+% ylabel('CR [dB]')
+% fprintf('CR: %f\n',[CRis(1)]);
+
+% show reference patch
+% figure 
+% imagesc(R21)
+
+end
\ No newline at end of file

Addons/CropZeros.m

+function OutMat = CropZeros(InMat)
+    % CropZeros removes border rows and columns if they have too many zeros.
+    %
+    % Input:
+    %   InMat - The input matrix.
+    %
+    % Output:
+    %   OutMat - The cropped matrix.
+
+    % Define the threshold for zero ratio
+    zeroThreshold = 0.95;
+
+    % Get matrix dimensions
+    [rows, cols] = size(InMat);
+
+    % Initialize start and end indices for rows and columns
+    rowStart = 1;
+    rowEnd = rows;
+    colStart = 1;
+    colEnd = cols;
+
+    % Check the top rows for zero ratio
+    while rowStart <= rowEnd && nnz(InMat(rowStart, :) == 0) / cols > zeroThreshold
+        rowStart = rowStart + 1;
+    end
+
+    % Check the bottom rows for zero ratio
+    while rowStart <= rowEnd && nnz(InMat(rowEnd, :) == 0) / cols > zeroThreshold
+        rowEnd = rowEnd - 1;
+    end
+
+    % Check the left columns for zero ratio
+    while colStart <= colEnd && nnz(InMat(:, colStart) == 0) / rows > zeroThreshold
+        colStart = colStart + 1;
+    end
+
+    % Check the right columns for zero ratio
+    while colStart <= colEnd && nnz(InMat(:, colEnd) == 0) / rows > zeroThreshold
+        colEnd = colEnd - 1;
+    end
+
+    % Crop the matrix based on the computed indices
+    OutMat = InMat(rowStart:rowEnd, colStart:colEnd);
+end

Addons/IsPatchValid.m

+function isValid = IsPatchValid(InMat)
+    % IsPatchValid checks if all quadrants of InMat contain at least one non-zero element.
+    %
+    % Input:
+    %   InMat - The input matrix.
+    %
+    % Output:
+    %   isValid - Returns true if all quadrants have at least one non-zero element; false otherwise.
+
+    % Get the dimensions of the matrix
+    [rows, cols] = size(InMat);
+
+    % Find the midpoint indices for rows and columns
+    midRow = ceil(rows / 2);
+    midCol = ceil(cols / 2);
+
+    % Define the four quadrants
+    Q1 = InMat(1:midRow, 1:midCol);           % Top-left quadrant
+    Q2 = InMat(1:midRow, midCol+1:end);       % Top-right quadrant
+    Q3 = InMat(midRow+1:end, 1:midCol);       % Bottom-left quadrant
+    Q4 = InMat(midRow+1:end, midCol+1:end);   % Bottom-right quadrant
+
+    % Check if each quadrant contains at least one non-zero element
+    isValid = any(Q1(:)) && any(Q2(:)) && any(Q3(:)) && any(Q4(:));
+end

Addons/WriteTif.m

+function WriteTif(MatIn,OutPutColorMap,tif_filename,varargin)
+%% function WriteTif(MatIn,OutPutColorMap,tif_filename,varargin)
+% created by Arthur Chavignon 18/12/19
+% example WriteTif(MatOutConv(:,:,1:10:end),hot(256),'tiftest.tif')
+
+if ~strcmp(tif_filename(end-3:end),'.tif')
+    error('OutPutName should end by .tif ')
+end
+
+tmp = strcmpi(varargin,'overwrite');
+if any(tmp),ForceOverWirte= varargin{find(tmp)+1};else, ForceOverWirte=0;end
+
+if exist(tif_filename,'file')==2
+    if ForceOverWirte==1
+        choice = 'Overwrite';
+    else
+        % Ask the user whether to overwrite the existing file:
+        choice = questdlg(['The file  ',tif_filename,' already exists. Overwrite it?'], ...
+            'The file already exists.','Overwrite','Cancel','Cancel');
+    end
+    % Overwriting basically means deleting and starting from scratch:
+    
+    if strcmp(choice,'Overwrite')
+        delete(tif_filename)
+    else
+        clear tif_filename firstframe DelayTime DitherOption LoopCount frame
+        error('The tiffing has been canceled.')
+    end
+end
+
+c0 = [0 max(MatIn(:))];
+tmp = strcmpi(varargin,'caxis');
+if any(tmp)
+    InputCaxis = varargin{find(tmp)+1};
+    if numel(InputCaxis)==1
+        c0 = c0.*InputCaxis;
+    elseif numel(InputCaxis)==2
+        c0 = InputCaxis;
+    end
+end
+
+%% Convert into digits
+if ~isempty(OutPutColorMap)
+    Ndigit = max(size(OutPutColorMap));
+    
+    MatInSat = MatIn;
+    MatInSat(MatIn>c0(2)) = c0(2);
+    MatInSat(MatInSat<c0(1)) = c0(1);
+    MatInNorm = [MatInSat-c0(1)]/diff(c0);
+    
+    MatInNorm = round(MatInNorm*Ndigit);
+    MatInNorm(MatInNorm>Ndigit)=Ndigit;
+else
+    % if MatIn is already RGB
+    MatInNorm = MatIn;
+end
+% Write RGB data
+
+for ii=1:size(MatInNorm,3)
+    if ~isempty(OutPutColorMap)
+        RGB = ind2rgb(MatInNorm(:,:,ii),OutPutColorMap); % convert to RGB
+    else
+        RGB = squeeze(MatInNorm(:,:,ii,:));
+    end
+    imwrite(RGB, tif_filename,'WriteMode','append','Compression','lzw','Description','ArthurChavignon');
+end
+
+tmp = strcmpi(varargin,'VoxelSizeMm');
+if any(tmp),VoxelSizeMm= varargin{find(tmp)+1};
+    t_tif = Tiff(tif_filename,'r+');
+    VoxSize_cm = VoxelSizeMm*1e-3;
+    setTag(t_tif,'XResolution',1/VoxSize_cm(1))
+    setTag(t_tif,'YResolution',1/VoxSize_cm(2))
+    setTag(t_tif,'ResolutionUnit',Tiff.ResolutionUnit.Centimeter)
+    close(t_tif);
+end
+
+return
\ No newline at end of file

Addons/blockproc3.m

+function im2 = blockproc3(im, blksz, fun, border, numworkers)
+% BLOCKPROC3  Block processing for a 3D image
+%
+% Sometimes images are too large to be processed in one block. Matlab
+% provides function blockproc(), that works on 2D images only. This
+% blockproc3() function, on the other hand, works also with 3D images.
+%
+% IM2 = blockproc3(IM, BLKSZ, FUN)
+%
+%   IM is a 2D or 3D-array with the input grayscale image.
+%
+%   IM2 is the ouput image, processed by blocks, with the same size as IM.
+%
+%   BLKSZ is a 2- or 3-vector with the size of the blocks the image will be
+%   split into. The last block in every row, column and slice will be
+%   smaller, if the image cannot accommodate it.
+%
+%   FUN is a function handle. This is the processing applied to every
+%   block.
+%
+% IM2 = blockproc3(IM, BLKSZ, FUN, BORDER)
+%
+%  BORDER is a 2- or 3-vector with the size of the border around each
+%  block. Basically, each block size will be grown by this amount, the
+%  block processed, and then the border cut off from the solution. This is
+%  useful to avoid "seams" in IM2.
+%
+% Example:
+%
+%       fun = @(x) deconvblind(x, ones(20, 20, 10));
+%       im2 = blockproc3d(im, [256 256 128], fun, [30 30 10]);
+%
+% IM2 = blockproc3(..., NUMWORKERS)
+%
+%   NUMWORKERS is an integer to activate parallel processing. If the
+%   Parallel Computing Toolbox is available and NUMWORKERS>1, each image
+%   block will be processed as a parallel job, using up to NUMWORKERS cores
+%   at the same time. By default, NUMWORKERS=1, and parallel processing is
+%   disabled.
+%
+%   Note that in parallel model, FUN gets wrapped in parallel jobs. Thus,
+%   execution can be interrupted by pressing CTRL+C even when FUN is a C++
+%   MEX file.
+%
+%   Parallel processing works by first splitting the image into blocks with
+%   borders, assigning processing of each block to a separate job, then
+%   processing each block, and finally trimming the borders and
+%   reassembling the results into the output image.
+%
+%   This will increase the amount of necessary memory, but this function is
+%   relatively efficient in the sense that only blocks being processed take
+%   up extra memory (blocks waiting in the job queue don't). 
+%
+%   The speed improvement may not be as good as expected (e.g. 4 processor
+%   will make the processing only 2x faster), as Matlab function can be
+%   already more or less optimised to make use of several processors.
+%
+%   To use this option, it is important to NOT have created a pool of
+%   workers in Matlab with "matlabpool open", as this would block all local
+%   resources, and not leave any cores free for processing the image
+%   blocks. A simple example:
+%
+%      % median filtering using a 19x19x19 voxel neighbourhood, processing
+%      % the image by blocks, using parallel computations
+%      matlabpool open
+%      sz = [19 19 19];
+%      fun = @(x) medfilt3(x, sz);
+%      im2 = blockproc3(im, [128 128 64], fun, (sz+1)/2, true);
+%      matlabpool close
+%
+%
+% See also: blockproc, scimat_blockproc3.
+
+% Author: Ramon Casero <rcasero@gmail.com>
+% Copyright © 2011 University of Oxford
+% Version: 0.4.0
+% 
+% University of Oxford means the Chancellor, Masters and Scholars of
+% the University of Oxford, having an administrative office at
+% Wellington Square, Oxford OX1 2JD, UK. 
+%
+% This file is part of Gerardus.
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details. The offer of this
+% program under the terms of the License is subject to the License
+% being interpreted in accordance with English Law and subject to any
+% action against the University of Oxford being under the jurisdiction
+% of the English Courts.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+% check arguments
+narginchk(3, 5);
+nargoutchk(0, 1);
+
+% defaults
+if isempty(blksz)
+    blksz = size(im);
+end
+if (nargin < 4) || isempty(border)
+    border = [0 0 0];
+end
+if (nargin < 5) || isempty(numworkers)
+    numworkers = 1;
+end
+
+% for convenience, we need the size vector to have 3 components, even for a
+% 2D image
+if (length(blksz) < 3)
+    blksz(3) = 1;
+end
+
+% image size
+imsz = size(im);
+
+% block limits without the extra borders...
+
+% ... starting points
+r0 = 1:blksz(1):imsz(1);
+c0 = 1:blksz(2):imsz(2);
+s0 = 1:blksz(3):imsz(3);
+
+% ... end points
+rx = r0 + blksz(1) - 1;
+cx = c0 + blksz(2) - 1;
+sx = s0 + blksz(3) - 1;
+rx = min(rx, imsz(1));
+cx = min(cx, imsz(2));
+sx = min(sx, imsz(3));
+
+% block limits with the extra borders...
+
+% ... starting points
+br0 = max(r0 - border(1), 1);
+bc0 = max(c0 - border(2), 1);
+bs0 = max(s0 - border(3), 1);
+
+% ... end points
+brx = min(rx + border(1), imsz(1));
+bcx = min(cx + border(2), imsz(2));
+bsx = min(sx + border(3), imsz(3));
+
+% number of blocks in each dimension
+NR = length(r0);
+NC = length(c0);
+NS = length(s0);
+
+% init output
+im2 = im;
+
+
+if (numworkers > 1) % parallel processing
+    
+    % build a cluster from the local profile
+    c = parcluster('local');
+    
+    if (c.NumWorkers < numworkers)
+        warning(['Number of workers selected by user is too large. Only ' ...
+            num2str(c.NumWorkers) ' cores available'])
+    else
+        % save number of maximum number of workers to undo the change after
+        % the function exits
+        numworkers0 = c.NumWorkers;
+        
+        % change profile's maximum number of workers to the number selected
+        % by the user
+        c.NumWorkers = numworkers;
+        
+        % when the function exits for whatever reason, reset the number of
+        % available workers
+        numworkersCleanupTrigger = onCleanup(@() resetNumworkers(numworkers0));
+    end
+    
+    % number of blocks
+    numblocks = NR * NC * NS;
+    
+    % init cell array to contain the job descriptors and job cleanup triggers
+    job = cell(1, numblocks);
+    jobCleanupTrigger = cell(1, numblocks);
+    
+    % process all input blocks (loops are in inverted order, so that
+    % linear indices follow 1, 2, 3, 4...)
+    for B = 1:numblocks
+        
+        % get r, c, s indices for current block
+        [I, J, K] = ind2sub([NR, NC, NS], B);
+        
+        % create a new job for this block
+        job{B} = createJob(c);
+        
+        % link the job to cleanup, so that if the user interrupts with
+        % CTRL+C, the job is cancelled and doesn't keep consuming resources
+        %
+        % note that we need to create a vector of onCleanup objects,
+        % because onCleanup gets called when the object is destroyed. and
+        % this will happen if we overwrite the same variable
+        jobCleanupTrigger{B} = onCleanup(@() cleanup(job{B}));
+ 
+        % add processing of current block as a task
+        createTask(job{B}, fun, 1, ...
+            {im(br0(I):brx(I), ...
+            bc0(J):bcx(J), ...
+            bs0(K):bsx(K) ...
+            )});
+        
+        % run job
+        submit(job{B});
+            
+    end
+    
+    % wait on all the jobs to finish, and extract the outputs
+    for B = 1:numblocks
+        
+        % get r, c, s indices for current block
+        [I, J, K] = ind2sub([NR, NC, NS], B);
+        
+        % wait until job finishes
+        wait(job{B});
+        
+        % extract the output of this job
+        aux = fetchOutputs(job{B});
+        
+        % assign result to output removing the borders
+        im2(...
+            r0(I):rx(I), ...
+            c0(J):cx(J), ...
+            s0(K):sx(K) ...
+            ) ...
+            = aux{1}(...
+            r0(I)-br0(I)+1:rx(I)-br0(I)+1, ...
+            c0(J)-bc0(J)+1:cx(J)-bc0(J)+1, ...
+            s0(K)-bs0(K)+1:sx(K)-bs0(K)+1 ...
+            );
+
+    end
+    
+else % single processor (we save memory by not creating a cell vector with all the blocks)
+
+    % iterate all image blocks
+    for I = 1:NR
+        for J = 1:NC
+            for K = 1:NS
+
+                % process current image block
+                aux = fun(im(...
+                    br0(I):brx(I), ...
+                    bc0(J):bcx(J), ...
+                    bs0(K):bsx(K) ...
+                    ));
+                
+                
+                % assign result to output removing the borders
+                im2(...
+                    r0(I):rx(I), ...
+                    c0(J):cx(J), ...
+                    s0(K):sx(K) ...
+                    ) ...
+                    = aux(...
+                    (r0(I)-br0(I)+1):(end-(brx(I)-rx(I))), ...
+                    (c0(J)-bc0(J)+1):(end-(bcx(J)-cx(J))), ...
+                    (s0(K)-bs0(K)+1):(end-(bsx(K)-sx(K))) ...
+                    );
+                
+            end
+        end
+    end
+
+end
+
+end
+
+% cleanup(): function called when we exit, in charge of cancelling one of
+% the current jobs. This can be due to two reasons: 1) the task finished
+% successfully or 2) the user interrupted with CTRL+C
+function cleanup(job)
+
+switch job.State
+    case {'running', 'queued'}
+        % cancel job so that it doesn't keep taking up resources
+        job.cancel;
+    case 'finished'
+        % DEBUG:
+%         disp('my_cleanup called because job finished');
+end
+
+end % function cleanup()
+
+% resetNumworkers(): function called when we exit, in charge of leaving the
+% cluster profile with the same number of maximum workers it had when we
+% entered the function
+function resetNumworkers(numworkers0)
+
+    % get cluster from the local profile
+    c = parcluster('local');
+    
+    % reset maximum number of workers
+    c.NumWorkers = numworkers0;
+
+end