Permalink
Browse files

- initial commit

  • Loading branch information...
mstorath committed Apr 20, 2017
0 parents commit 452b1cecd4998c3725ace36ebe3e452f8464833c
Showing with 4,431 additions and 0 deletions.
  1. +30 −0 Auxiliary/Operators/convop.m
  2. +66 −0 Auxiliary/Operators/linop.m
  3. +25 −0 Auxiliary/Operators/radonop.m
  4. +46 −0 Auxiliary/convkernel.m
  5. +6 −0 Auxiliary/countJumps.m
  6. +66 −0 Auxiliary/energyL2Potts.m
  7. +13 −0 Auxiliary/intmatrix.m
  8. +21 −0 Auxiliary/medianw.m
  9. +7 −0 Auxiliary/plimage2double.m
  10. +10 −0 Auxiliary/plpsnr.m
  11. +13 −0 Auxiliary/psnr.m
  12. +11 −0 Auxiliary/randidx.m
  13. +11 −0 Auxiliary/randl.m
  14. +20 −0 Auxiliary/samplesToWeights.m
  15. +40 −0 Auxiliary/setPLJavaPath.m
  16. +30 −0 Auxiliary/showPotts.m
  17. +34 −0 Auxiliary/showSparse.m
  18. +46 −0 Auxiliary/showdbl.m
  19. +8 −0 Auxiliary/softThreshold.m
  20. +21 −0 Auxiliary/spconvmatrix.m
  21. +12 −0 Auxiliary/spdiffmatrix.m
  22. BIN Data/church.jpg
  23. BIN Data/colors.png
  24. BIN Data/desert.jpg
  25. +11 −0 Data/imageSources.txt
  26. +151 −0 Data/loadPcwConst.m
  27. +102 −0 Data/loadSparse.m
  28. BIN Data/lookmickey.jpg
  29. +33 −0 Demos/1D/demoL1iPotts_DecImp.m
  30. +32 −0 Demos/1D/demoL1iPotts_DecLap.m
  31. +30 −0 Demos/1D/demoL1iSpars_DecImp.m
  32. +27 −0 Demos/1D/demoL1iSpars_DecLap.m
  33. +26 −0 Demos/1D/demoL2iPotts_Deconv.m
  34. +29 −0 Demos/1D/demoL2iPotts_Fourier.m
  35. +34 −0 Demos/1D/demoL2iPotts_Linop.m
  36. +27 −0 Demos/1D/demoL2iSpars_Deconv.m
  37. +19 −0 Demos/1D/demoPotts_GaussianNoise.m
  38. +21 −0 Demos/1D/demoPotts_ImpulsiveNoise.m
  39. +19 −0 Demos/1D/demoPotts_LaplacianNoise.m
  40. +25 −0 Demos/2D/demoPotts2DColorDenoising.m
  41. +33 −0 Demos/2D/demoPotts2DColorInpainting.m
  42. +32 −0 Demos/2D/demoPotts2DColorMissing.m
  43. +22 −0 Demos/2D/demoPotts2DColorSegmentation.m
  44. +39 −0 Demos/2D/demoPotts2DRadon.m
  45. BIN Demos/2D/results/plFigPottsRec7Angles.fig
  46. +8 −0 Demos/pottsLabDemo.m
  47. BIN Java/bin/pottslab/IndexedLinkedHistogram$HistNode.class
  48. BIN Java/bin/pottslab/IndexedLinkedHistogram.class
  49. BIN Java/bin/pottslab/IndexedLinkedHistogramUnweighted$HistNode.class
  50. BIN Java/bin/pottslab/IndexedLinkedHistogramUnweighted.class
  51. BIN Java/bin/pottslab/JavaTools.class
  52. BIN Java/bin/pottslab/L2Potts.class
  53. BIN Java/bin/pottslab/PLImage$1.class
  54. BIN Java/bin/pottslab/PLImage.class
  55. BIN Java/bin/pottslab/PLProcessor.class
  56. BIN Java/bin/pottslab/PLVector.class
  57. +342 −0 Java/src/pottslab/IndexedLinkedHistogram.java
  58. +240 −0 Java/src/pottslab/IndexedLinkedHistogramUnweighted.java
  59. +309 −0 Java/src/pottslab/JavaTools.java
  60. +136 −0 Java/src/pottslab/L2Potts.java
  61. +142 −0 Java/src/pottslab/PLImage.java
  62. +306 −0 Java/src/pottslab/PLProcessor.java
  63. +258 −0 Java/src/pottslab/PLVector.java
  64. +116 −0 Potts/PottsCore/findBestPartition.m
  65. +147 −0 Potts/PottsCore/iPottsADMM.m
  66. +56 −0 Potts/PottsCore/reconstructionFromPartition.m
  67. +80 −0 Potts/minL1Potts.m
  68. +22 −0 Potts/minL1iPotts.m
  69. +55 −0 Potts/minL2Potts.m
  70. +22 −0 Potts/minL2iPotts.m
  71. +203 −0 Potts2D/Core/iPotts2DADMM.m
  72. +68 −0 Potts2D/minL2Potts2DADMM.m
  73. +5 −0 Potts2D/minL2iPotts2DADMM.m
  74. +71 −0 Readme.txt
  75. +138 −0 Sparsity/SparsityCore/iSparsADMM.m
  76. +54 −0 Sparsity/SparsityCore/iSparsByPottsADMM.m
  77. +14 −0 Sparsity/SparsityCore/minSpars.m
  78. +11 −0 Sparsity/minL1Spars.m
  79. +25 −0 Sparsity/minL1iSpars.m
  80. +11 −0 Sparsity/minL2Spars.m
  81. +36 −0 Sparsity/minL2iSpars.m
  82. +117 −0 Tikhonov/minL1Tikhonov.m
  83. +98 −0 Tikhonov/minL2Tikhonov.m
  84. +57 −0 Tikhonov/minL2TikhonovFBP.m
  85. +36 −0 installPottslab.m
@@ -0,0 +1,30 @@
+classdef convop < double
+%convop A class for circular convolution by FFT
+
+% written by M. Storath
+% $Date: 2014-05-07 11:42:11 +0200 (Mi, 07. Mai 2014) $ $Revision: 89 $
+
+
+ methods
+ % constructor
+ function obj = convop(fourier)
+ obj = obj@double(fourier);
+ end
+
+ % ctranspose (')
+ function C = ctranspose(A)
+ C = convop(conj(A));
+ end
+
+ % mtimes (*)
+ function C = mtimes( A, B )
+ if isa(B, 'convop')
+ C = convop( A .* B );
+ else
+ C = ifftn( A .* fftn(B) );
+ end
+ end
+
+ end
+end
+
@@ -0,0 +1,66 @@
+classdef linop
+%linop A class for linear operators
+
+% written by M. Storath
+% $Date: 2015-10-15 15:07:43 +0200 (Do, 15. Okt 2015) $ $Revision: 132 $
+
+
+ properties
+ % function handle for evaluation A * x
+ eval;
+ % function handle for evaluation A' * x
+ ctrans;
+ % A'A
+ normalOp;
+ % true if A'A positive definite
+ posdef;
+ % true if A'A positive definite
+ lseSolver;
+ end
+
+ methods
+ % constructor
+ function A = linop(eval, ctrans, varargin)
+ % if it is a matrix
+ if ~isa(eval, 'function_handle')
+ A.eval = @(x) eval * x;
+ A.ctrans = @(x) eval' * x;
+ M = eval' * eval;
+ A.normalOp = @(x) M * x;
+ else
+ A.eval = eval;
+ A.ctrans = ctrans;
+ A.normalOp = @(x) A.ctrans(A.eval(x));
+ end
+ ip = inputParser;
+ addParamValue(ip, 'posdef', false);
+ addParamValue(ip, 'lseSolver', []);
+ parse(ip, varargin{:});
+ par = ip.Results;
+ A.posdef = par.posdef;
+ A.lseSolver = par.lseSolver;
+ end
+
+ % ctranspose (')
+ function C = ctranspose(A)
+ C = linop(A.ctrans, A.eval);
+ end
+
+ % mtimes (*)
+ function C = mtimes( A, B )
+ if isa(B, 'linop')
+ C = linop( @(x) A.eval(B.eval(x)), @(x) B.ctrans(A.ctrans(x)));
+ else
+ C = A.eval(B);
+ end
+ end
+
+ % size
+ function s = size(A)
+ %warning('Size is deprecated for class linop');
+ s = [1 1];
+ end
+
+ end
+end
+
@@ -0,0 +1,25 @@
+classdef radonop < linop
+ %radonop A class for Radon transform
+
+ % written by M. Storath
+ % $Date: 2014-09-02 14:56:23 +0200 (Di, 02. Sep 2014) $ $Revision: 105 $
+
+ properties
+ % solution method for prox
+ useFBP;
+ end
+
+ methods
+ % constructor
+ function A = radonop(theta, imgSize)
+
+ eval = @(x) radon(x, theta);
+ ctrans = @(x) iradon(x, theta, 'linear', 'None', 1, imgSize);
+
+ A = A@linop(eval, ctrans);
+ A.posdef = true;
+ A.useFBP = false;
+ end
+ end
+end
+
@@ -0,0 +1,46 @@
+function h = convkernel( type, n, scale )
+%CONVKERNEL Creates a convolution kernel
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+if not(exist('scale', 'var'))
+ scale = 1;
+end
+
+% set coordinates
+d = n/2;
+x = linspace(-d, d, n)/scale;
+
+switch type
+ case 'box'
+ h = ones(n, 1);
+
+ case 'hat'
+ h(x >= 0) = 1-x(x >= 0);
+ h(x < 0) = 1 + x(x < 0);
+ h(abs(x) > 1) = 0;
+
+ case 'gaussian'
+ h = exp(-0.5 * x.^2);
+
+ case 'doubleexp'
+ h = exp(-0.5 * abs(x));
+
+ case 'mollifier'
+ h = exp(-1./(1 - abs(x).^2));
+ h(abs(x) > 1) = 0;
+
+ case 'ramp'
+ h(x >= 0) = 1-x(x >= 0);
+ h(abs(x) > 1) = 0;
+
+ otherwise
+ error('This option does not exist.')
+end
+
+% normalize
+h = h(:)./sum(h);
+
+end
+
@@ -0,0 +1,6 @@
+function j = countJumps( f )
+%COUNTJUMPS Counts the jumps of a vector f
+j = sum( diff(f) ~= 0);
+
+end
+
@@ -0,0 +1,66 @@
+function [energy, jumpPenalty, dataError] = energyL2Potts( u, f, gamma, A, isotropic )
+%energyL2Potts Computes the energy of the L2 Potts functional
+
+if exist('A', 'var') && not(isempty(A))
+ Au = A * u;
+ dataError = sum((Au(:) - f(:)).^2);
+else
+ dataError = sum((u(:) - f(:)).^2);
+end
+
+
+if isvector(u)
+ % vectors
+ jumpPenalty = sum(diff(u(:)) ~= 0);
+else
+ % matrices
+ nJumpComp = 0; % jumps in compass directions
+ nJumpDiag = 0; % jumps in diagonal directions
+
+ % count jumps
+ for i = 1:size(u, 1)
+ for j = 1:size(u, 2)-1
+ if ~all((u(i,j,:) == u(i,j+1,:)))
+ nJumpComp = nJumpComp + 1;
+ end
+ end
+ end
+ for i = 1:size(u, 1)-1
+ for j = 1:size(u, 2)
+ if ~all((u(i,j,:) == u(i+1,j,:)))
+ nJumpComp = nJumpComp + 1;
+ end
+ end
+ end
+ for i = 1:size(u, 1)-1
+ for j = 1:size(u, 2)-1
+ if ~all((u(i,j,:) == u(i+1,j+1,:)))
+ nJumpDiag = nJumpDiag + 1;
+ end
+ end
+ end
+ for i = 1:size(u, 1)-1
+ for j = 2:size(u, 2)
+ if ~all((u(i,j,:) == u(i+1,j-1,:)))
+ nJumpDiag = nJumpDiag + 1;
+ end
+ end
+ end
+
+ % set weights (isotropic by default)
+ if ~exist('isotropic', 'var') || isotropic
+ omega1 = sqrt(2) - 1;
+ omega2 = 1 - sqrt(2)/2;
+ else
+ omega1 = 1;
+ omega2 = 0;
+ end
+
+ % compute energy
+ jumpPenalty = (omega1 * nJumpComp + omega2* nJumpDiag);
+
+end
+
+energy = gamma * jumpPenalty + dataError;
+
+end
@@ -0,0 +1,13 @@
+function A = intmatrix( n )
+%INTMATRIX Creates a matrix that acts as integration
+% (cumulative sum) on the input vector
+%
+% Au = cumsum(u)
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+A = tril(ones(n));
+
+end
+
@@ -0,0 +1,21 @@
+function mu = medianw( data, weight )
+% medianw Computes a weighted median of the data
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+if numel(data) == 1
+ mu = data;
+else
+ [dataSorted, idx] = sort(data);
+ weightSorted = weight(idx);
+ weightSum = sum(weight);
+ cumWeight = cumsum(weightSorted)/weightSum;
+ % find first median index
+ idx = find(cumWeight < 0.5, 1, 'last' ) + 1;
+ if isempty(idx)
+ idx = 1;
+ end
+ mu = dataSorted(idx);
+end
+
@@ -0,0 +1,7 @@
+function img = plimage2double( plimg )
+%PLIMAGE2DOUBLE Summary of this function goes here
+% Detailed explanation goes here
+img = reshape(plimg.toDouble(), [plimg.mRow, plimg.mCol, plimg.mLen]);
+
+end
+
@@ -0,0 +1,10 @@
+function r = plpsnr( groundTruth, signal )
+%plPSNR Computes the peak signal to noise ratio
+
+% written by M. Storath
+% $Date: 2014-05-07 14:23:10 +0200 (Mi, 07. Mai 2014) $ $Revision: 90 $
+
+mse = mean(abs(signal(:) - groundTruth(:)).^2);
+r = 10 * log10(max(groundTruth(:))^2 / mse);
+
+end
@@ -0,0 +1,13 @@
+function r = psnr( groundTruth, signal )
+%PSNR Computes the peak signal to noise ratio
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+warning('The function psnr conflicts with the new function psnr from the image processing toolbox from Matlab v2014a. Use plpsnr instead.')
+
+mse = mean(abs(signal(:) - groundTruth(:)).^2);
+r = 10 * log10(max(groundTruth(:))^2 / mse);
+
+end
+
@@ -0,0 +1,11 @@
+function idx = randidx( siz, fraction )
+%randidx Selects randomly a fraction of indices for a matrix of size siz
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+n = prod(siz);
+idx = randperm(n) ;
+idx = idx(1:round(fraction * n));
+end
+
@@ -0,0 +1,11 @@
+function y = randl( varargin )
+%randl Random normal Laplacian noise, mean = 0, variance sigma^2 = 1
+% pdf: $\frac{1}{\sqrt{2}} e^{- \sqrt{2} |x|}$
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+x = rand( varargin{:} ) - 0.5;
+y = -sign(x) .* log(1 - 2 * abs(x)) / sqrt(2);
+end
+
@@ -0,0 +1,20 @@
+function weights = samplesToWeights( samples )
+%samplesToWeights Computes weights from a set of sample points
+
+% written by M. Storath
+% $Date: 2012/10/29 01:19:08 $ $Revision: 0.1 $
+
+%init
+weights = zeros(size(samples));
+
+% first and last weight
+weights(1) = samples(2) - samples(1);
+weights(end) = samples(end) - samples(end-1);
+
+% central weights
+subs = 2:numel(samples)-1;
+weights(subs) = 0.5 * ( abs(samples(subs) - samples(subs-1)) ...
+ + abs(samples(subs) - samples(subs+1)) );
+
+end
+
Oops, something went wrong.

0 comments on commit 452b1ce

Please sign in to comment.