diff --git a/Addons/ContrastRatio.m b/Addons/ContrastRatio.m
new file mode 100644
index 0000000000000000000000000000000000000000..b92c1220e7c27d6d0b18191881c8456a0772461b
--- /dev/null
+++ b/Addons/ContrastRatio.m
@@ -0,0 +1,59 @@
+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
diff --git a/Addons/CropZeros.m b/Addons/CropZeros.m
new file mode 100644
index 0000000000000000000000000000000000000000..ff833d022b76c0c2a9902d1ef6aec7aae5598539
--- /dev/null
+++ b/Addons/CropZeros.m
@@ -0,0 +1,44 @@
+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
diff --git a/Addons/IsPatchValid.m b/Addons/IsPatchValid.m
new file mode 100644
index 0000000000000000000000000000000000000000..63858ca401e4966941442b1136c7125e72ad4405
--- /dev/null
+++ b/Addons/IsPatchValid.m
@@ -0,0 +1,25 @@
+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
diff --git a/Addons/WriteTif.m b/Addons/WriteTif.m
new file mode 100644
index 0000000000000000000000000000000000000000..1de411a1481479a3f3bfa93e89817a69fd074ca7
--- /dev/null
+++ b/Addons/WriteTif.m
@@ -0,0 +1,78 @@
+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
diff --git a/Addons/blockproc3.m b/Addons/blockproc3.m
new file mode 100644
index 0000000000000000000000000000000000000000..b48f99233a9e6e410cc5b91138cf98f97c14c665
--- /dev/null
+++ b/Addons/blockproc3.m
@@ -0,0 +1,308 @@
+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