From f618466c25d43f3bae9e40920273bf77de1e1149 Mon Sep 17 00:00:00 2001 From: leochanj105 Date: Mon, 19 Oct 2020 23:09:30 -0400 Subject: initial sd-vbs initial sd-vbs add sd-vbs sd-vbs --- SD-VBS/common/toolbox/MultiNcut/MNcut.m | 93 ++ SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m | 34 + SD-VBS/common/toolbox/MultiNcut/README.tex | 9 + SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c | 405 ++++++++ .../common/toolbox/MultiNcut/a_times_b_cmplx.dll | Bin 0 -> 7114 bytes .../toolbox/MultiNcut/a_times_b_cmplx.mexa64 | Bin 0 -> 10646 bytes .../toolbox/MultiNcut/a_times_b_cmplx.mexglx | Bin 0 -> 7896 bytes .../toolbox/MultiNcut/a_times_b_cmplx.mexmac | Bin 0 -> 13096 bytes SD-VBS/common/toolbox/MultiNcut/affinityic.c | 187 ++++ SD-VBS/common/toolbox/MultiNcut/affinityic.dll | Bin 0 -> 7680 bytes SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64 | Bin 0 -> 9755 bytes SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx | Bin 0 -> 7775 bytes SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m | 48 + SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c | 198 ++++ SD-VBS/common/toolbox/MultiNcut/cimgnbmap.dll | Bin 0 -> 7168 bytes SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexa64 | Bin 0 -> 9598 bytes SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexglx | Bin 0 -> 7824 bytes SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.c | 197 ++++ .../common/toolbox/MultiNcut/cimgnbmap_cross.dll | Bin 0 -> 7168 bytes .../toolbox/MultiNcut/cimgnbmap_cross.mexa64 | Bin 0 -> 9604 bytes .../toolbox/MultiNcut/cimgnbmap_cross.mexglx | Bin 0 -> 7830 bytes SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.c | 294 ++++++ SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.dll | Bin 0 -> 7680 bytes .../common/toolbox/MultiNcut/cimgnbmap_star.mexa64 | Bin 0 -> 9747 bytes .../common/toolbox/MultiNcut/cimgnbmap_star.mexglx | Bin 0 -> 7829 bytes SD-VBS/common/toolbox/MultiNcut/compileAll.m | 10 + SD-VBS/common/toolbox/MultiNcut/computeMultiW.m | 245 +++++ SD-VBS/common/toolbox/MultiNcut/discretisation.m | 49 + .../MultiNcut/discretisationEigenVectorData.m | 12 + SD-VBS/common/toolbox/MultiNcut/doog1.m | 32 + SD-VBS/common/toolbox/MultiNcut/doog2.m | 38 + SD-VBS/common/toolbox/MultiNcut/eigSolve.m | 5 + SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m | 29 + SD-VBS/common/toolbox/MultiNcut/gaussian.m | 31 + .../toolbox/MultiNcut/make_filterbank_even2.m | 45 + .../toolbox/MultiNcut/make_filterbank_odd2.m | 46 + .../common/toolbox/MultiNcut/mex_projection_QR.c | 82 ++ .../common/toolbox/MultiNcut/mex_projection_QR.dll | Bin 0 -> 10240 bytes .../toolbox/MultiNcut/mex_projection_QR.mexa64 | Bin 0 -> 13067 bytes .../toolbox/MultiNcut/mex_projection_QR.mexglx | Bin 0 -> 10165 bytes .../MultiNcut/mex_projection_QR_symmetric.c | 83 ++ .../MultiNcut/mex_projection_QR_symmetric.dll | Bin 0 -> 10240 bytes .../MultiNcut/mex_projection_QR_symmetric.mexa64 | Bin 0 -> 13093 bytes .../MultiNcut/mex_projection_QR_symmetric.mexglx | Bin 0 -> 10175 bytes .../MultiNcut/mex_w_times_x_symmetric.mexglx | Bin 0 -> 8713 bytes .../MultiNcut/mex_w_times_x_symmetric.mexmac | Bin 0 -> 13396 bytes SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c | 216 ++++ .../common/toolbox/MultiNcut/multiAffinityic.dll | Bin 0 -> 8192 bytes .../toolbox/MultiNcut/multiAffinityic.mexa64 | Bin 0 -> 10138 bytes .../toolbox/MultiNcut/multiAffinityic.mexglx | Bin 0 -> 8134 bytes .../toolbox/MultiNcut/multiIntensityFirstLayer.c | 126 +++ .../toolbox/MultiNcut/multiIntensityFirstLayer.dll | Bin 0 -> 7168 bytes .../MultiNcut/multiIntensityFirstLayer.mexa64 | Bin 0 -> 9345 bytes .../MultiNcut/multiIntensityFirstLayer.mexglx | Bin 0 -> 6805 bytes .../common/toolbox/MultiNcut/multiIntensityWppc.c | 158 +++ .../toolbox/MultiNcut/multiIntensityWppc.dll | Bin 0 -> 7168 bytes .../toolbox/MultiNcut/multiIntensityWppc.mexa64 | Bin 0 -> 9493 bytes .../toolbox/MultiNcut/multiIntensityWppc.mexglx | Bin 0 -> 7277 bytes SD-VBS/common/toolbox/MultiNcut/quadedgep2.m | 188 ++++ .../common/toolbox/MultiNcut/quickNcutHardBiais2.m | 187 ++++ SD-VBS/common/toolbox/MultiNcut/read_data.m | 13 + SD-VBS/common/toolbox/MultiNcut/readimage.m | 15 + SD-VBS/common/toolbox/MultiNcut/run_script.m | 60 ++ SD-VBS/common/toolbox/MultiNcut/showmask.m | 65 ++ SD-VBS/common/toolbox/MultiNcut/sparsifyc.c | 232 +++++ SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll | Bin 0 -> 8704 bytes SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64 | Bin 0 -> 10322 bytes SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx | Bin 0 -> 8296 bytes SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac | Bin 0 -> 9004 bytes SD-VBS/common/toolbox/MultiNcut/spmtimesd.c | 141 +++ SD-VBS/common/toolbox/MultiNcut/spmtimesd.dll | Bin 0 -> 7168 bytes SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexa64 | Bin 0 -> 9282 bytes SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexglx | Bin 0 -> 7280 bytes SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexmac | Bin 0 -> 8888 bytes SD-VBS/common/toolbox/MultiNcut/tim_eigs.m | 1084 ++++++++++++++++++++ 75 files changed, 4657 insertions(+) create mode 100755 SD-VBS/common/toolbox/MultiNcut/MNcut.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/README.tex create mode 100755 SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.mexmac create mode 100755 SD-VBS/common/toolbox/MultiNcut/affinityic.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/affinityic.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_cross.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/cimgnbmap_star.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/compileAll.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/computeMultiW.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/discretisation.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/doog1.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/doog2.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/eigSolve.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/gaussian.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/make_filterbank_even2.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/make_filterbank_odd2.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_projection_QR_symmetric.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/mex_w_times_x_symmetric.mexmac create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiAffinityic.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/quadedgep2.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/read_data.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/readimage.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/run_script.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/showmask.m create mode 100755 SD-VBS/common/toolbox/MultiNcut/sparsifyc.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac create mode 100755 SD-VBS/common/toolbox/MultiNcut/spmtimesd.c create mode 100755 SD-VBS/common/toolbox/MultiNcut/spmtimesd.dll create mode 100755 SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexa64 create mode 100755 SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexglx create mode 100755 SD-VBS/common/toolbox/MultiNcut/spmtimesd.mexmac create mode 100755 SD-VBS/common/toolbox/MultiNcut/tim_eigs.m (limited to 'SD-VBS/common/toolbox/MultiNcut') diff --git a/SD-VBS/common/toolbox/MultiNcut/MNcut.m b/SD-VBS/common/toolbox/MultiNcut/MNcut.m new file mode 100755 index 0000000..5486080 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/MNcut.m @@ -0,0 +1,93 @@ +function [NcutDiscretes,eigenVectors,eigenValues] = MNcut(I,nsegs); +% +% [NcutDiscrete,eigenVectors,eigenValues] = MNcut(I,nsegs); +% +% + +[nr,nc,nb] = size(I); + +max_image_size = max(nr,nc); + +% modified by song, 06/13/2005 +% test parameters +if (1) % original settings + if (max_image_size>120) & (max_image_size<=500), + % use 3 levels, + data.layers.number=3; + data.layers.dist=3; + data.layers.weight=[3000,4000,10000]; + data.W.scales=[1,2,3];%[1,2,3]; + data.W.radius=[2,3,7];%[2,3,7]; + elseif (max_image_size >500), + % use 4 levels, + data.layers.number=4; + data.layers.dist=3; + data.layers.weight=[3000,4000,10000,20000]; + data.W.scales=[1,2,3,3]; + data.W.radius=[2,3,4,6]; + elseif (max_image_size <=120) + data.layers.number=2; + data.layers.dist=3; + data.layers.weight=[3000,10000]; + data.W.scales=[1,2]; + data.W.radius=[2,6]; + end +else % test setting + if (max_image_size>200) & (max_image_size<=500), + % use 3 levels, + data.layers.number=3; + data.layers.dist=3; + data.layers.weight=[3000,4000,10000]; + data.W.scales=[1,2,3];%[1,2,3]; + data.W.radius=[2,3,7];%[2,3,7]; + elseif (max_image_size >500), + % use 4 levels, + data.layers.number=4; + data.layers.dist=3; + data.layers.weight=[3000,4000,10000,20000]; + data.W.scales=[1,2,3,3]; + data.W.radius=[2,3,4,6]; + elseif (max_image_size <=200) + data.layers.number=2; + data.layers.dist=3; + data.layers.weight=[3000,10000]; + data.W.scales=[1,2]; + data.W.radius=[2,4]; + end + +end; + + +data.W.edgeVariance=0.1; %0.1 +data.W.gridtype='square'; +data.W.sigmaI=0.12;%0.12 +data.W.sigmaX=1000; +data.W.mode='mixed'; +data.W.p=0; +data.W.q=0; + +%eigensolver +data.dataGraphCut.offset = 100;% 10; %valeur sur diagonale de W (mieux vaut 10 pour valeurs negatives de W) +data.dataGraphCut.maxiterations=50;% voir +data.dataGraphCut.eigsErrorTolerance=1e-2;%1e-6; +data.dataGraphCut.valeurMin=1e-6;%1e-5;% utilise pour tronquer des valeurs et sparsifier des matrices +data.dataGraphCut.verbose = 0; + +data.dataGraphCut.nbEigenValues=max(nsegs); + +disp('computeEdge'); +[multiWpp,ConstraintMat, Wind,data,emag,ephase]= computeMultiW (I,data); + +disp('Ncut'); +[eigenVectors,eigenValues]= eigSolve (multiWpp,ConstraintMat,data); + +%NcutDiscretes = zeros(nr,nc,length(nsegs)); +NcutDiscretes = zeros(nr,nc,(nsegs)); + +for j=1:length(nsegs), + nseg = nsegs(j); + [nr,nc,nb] = size(eigenVectors(:,:,1:nseg)); + [NcutDiscrete,evrotated] =discretisation(reshape(eigenVectors(:,:,1:nb),nr*nc,nb),nr,nc); + NcutDiscretes(:,:,j) = NcutDiscrete; +end + diff --git a/SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m b/SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m new file mode 100755 index 0000000..972a4eb --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/MNcutDemo.m @@ -0,0 +1,34 @@ +% MNcutDemo.m +% created by song, 06/13/2005 +% an exmaple of how to use and display MNcut + +num_segs = [20]; +imageSize = 800; + +img_filename = '/u/ikkjin/Benchmark/stitch/data/test/capitol/img1.jpg'; + +I=readimage(img_filename,imageSize); + +[SegLabel,eigenVectors,eigenValues]= MNcut(I,num_segs); + +for j=1:size(SegLabel,3), + [gx,gy] = gradient(SegLabel(:,:,j)); + bw = (abs(gx)>0.1) + (abs(gy) > 0.1); + + figure(1);clf; J1=showmask(double(I),bw); imagesc(J1);axis image; axis off; + set(gca, 'Position', [0 0 1 1]); + + % cm = sprintf('print -djpeg %s/file%.4d-%.2d.jpg',OutputDir,id,num_segs(j)); disp(cm);eval(cm); + + + % figure(10);imagesc(SegLabel(:,:,j));axis image; axis off; + % set(gca, 'Position', [0 0 1 1]); + % cm = sprintf('print -djpeg %s/Seg%.4d-%.2d.jpg',OutputDir,id,num_segs(j));disp(cm);eval(cm); + + % pause; +end + +% fname = files(id).name; +%cm = sprintf('save %s/SegLabl%.4d.mat I SegLabel fname',OutputDir,id); disp(cm); eval(cm); +%cm = sprintf('save %s/SegEig%.4d.mat eigenVectors eigenValues',OutputDir,id);disp(cm); eval(cm); + diff --git a/SD-VBS/common/toolbox/MultiNcut/README.tex b/SD-VBS/common/toolbox/MultiNcut/README.tex new file mode 100755 index 0000000..5970fb2 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/README.tex @@ -0,0 +1,9 @@ +1) You need to first compile the .c files,type + +>> compileAll('.'); + +2) the top level function is called MNcut.m + + + + diff --git a/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c new file mode 100755 index 0000000..25def92 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/a_times_b_cmplx.c @@ -0,0 +1,405 @@ +/*================================================================ +a_times_b_cmplx.c = used by a couple of mex functions +provide Matrix vector multiplications, +and solve triangular systems +(sparse matrix and full vector) + +CSC_CmplxVecMult_CAB_double, CSR_CmplxVecMult_CAB_double, +CSCsymm_CmplxVecMult_CAB_double added by Mirko Visontai (10/24/2003) + +*=================================================================*/ +# include "math.h" + + +/*C<-a*A*B+C*/ +void CSC_VecMult_CaABC_double( + const int m, const int k, const double alpha, + const double *val, const int *indx, + const int *pntrb, + const double *b, + double *c) +{ + int i,j,jb,je; + + for (i=0;i!=k;i++){ + jb = pntrb[i]; + je = pntrb[i+1]; + for (j=jb;j!=je;j++) + c[indx[j]] += alpha * b[i] * val[j]; + } +} + +/*C<-a*A'*B+C*/ +void CSR_VecMult_CaABC_double( + const int k, const int m, const double alpha, + const double *val, const int *indx, + const int *pntrb, + const double *b, + double *c) +{ + double t; + const double *pval; + int i,j,jb,je; + + pval = val; + for (i=0;i!=m;i++) { + t = 0; + jb = pntrb[i]; + je = pntrb[i+1]; + for (j=jb;j!=je;j++) + t += alpha * b[indx[j]] * (*pval++); + c[i] += t; + } +} + + +/*C<-A*b */ +void CSC_VecMult_CAB_double( + const int m, const int k, /*nb_rows, nb_columns*/ + const double *val, const int *indx, + const int *pntrb, + const double *b, + double *c + ) +{ + int i,j,jb,je; + double *pc=c; + for (i=0;i!=m;i++) *pc++ = 0; + + for (i=0;i!=k;i++){ + jb = pntrb[i]; + je = pntrb[i+1]; + for (j=jb;j!=je;j++) + c[indx[j]] += b[i] * val[j]; + } +} + +/*C<-A*b (complex)*/ +void CSC_CmplxVecMult_CAB_double( + const int m, const int k, + const double *valr, const double *vali, + const int *indx, + const int *pntrb, + const double *br, const double *bi, + double *cr, double *ci + ) +{ + int i,j,jb,je; + double *pcr=cr; + double *pci=ci; + for (i=0;i!=m;i++){ + *pcr++ = 0.0; + *pci++ = 0.0; + } + + for (i=0;i!=k;i++){ + jb = pntrb[i]; + je = pntrb[i+1]; + for (j=jb;j!=je;j++){ + cr[indx[j]] += (br[i] * valr[j]) - (bi[i] * vali[j]); + ci[indx[j]] += (br[i] * vali[j]) + (bi[i] * valr[j]); + } + } +} + +/*C<-A'*b + plus rapide que CSC_VecMult_CAB_double */ +void CSR_VecMult_CAB_double( + const int k, const int m, + const double *val, const int *indx, + const int *pntrb, + const double *b, + double *c + ) +{ + double t; + const double *pval; + double *pc=c; + int i,j,jb,je; + + for (i=0;i!=m;i++) *pc++ = 0; + + pval = val; + for (i=0;i!=m;i++) { + t = 0; + jb = pntrb[i]; + je = pntrb[i+1]; + for (j=jb;j!=je;j++) + t += b[indx[j]] * (*pval++); + c[i] += t; + } +} + +/*C<-A'*b (complex) + plus rapide que CSC_VecMult_CAB_double */ +void CSR_CmplxVecMult_CAB_double( + const int k, const int m, + const double *valr, const double *vali, + const int *indx, + const int *pntrb, + const double *br, const double *bi, + double *cr, double *ci + ) +{ + double tr, ti; + const double *pvalr; + const double *pvali; + double *pcr=cr; + double *pci=ci; + int i,j,jb,je; + + for (i=0;i!=m;i++){ + *pcr++ = 0.0; + *pci++ = 0.0; + } + + pvalr = valr; + pvali = vali; + for (i=0;i!=m;i++) { + tr = 0.0; + ti = 0.0; + jb = pntrb[i]; + je = pntrb[i+1]; + for (j=jb;j!=je;j++){ + tr += (br[indx[j]] * (*pvalr)) - (bi[indx[j]] * (*pvali)); + ti += (br[indx[j]] * (*pvali++)) + (bi[indx[j]] * (*pvalr++)); + } + cr[i] += tr; + ci[i] += ti; + } +} + + + +/* C<-A*b (A is symmetric) */ +void CSRsymm_VecMult_CAB_double( + const int k, const int m, + const double *val, const int *indx, + const int *pntrb, + const double *b, + double *c + ) +{ + const double *pval; + double *pc=c; + int i,j; + int jj; + int rpntrb, rpntre; + int index, nvals; + + + for (i=0;i!=m;i++) *pc++ = 0; + pval = val; + for (j=0;j!=k;j++){ + rpntrb = pntrb[j]; + rpntre = pntrb[j+1]; + for (jj=rpntrb;jj!=rpntre;jj++) { + index = indx[jj]; + if ( index == j ) { + c[j] += b[j] * (*pval++); + continue; + } + if ( index > j ) { + c[index] += b[j] * (*pval); + + c[j] += b[index] * (*pval++); + } + else { + pval++; + } + } + } +} + + +/* C<-A*b (A is symmetric and complex) */ +void CSRsymm_CmplxVecMult_CAB_double( + const int k, const int m, + const double *valr, const double *vali, + const int *indx, + const int *pntrb, + const double *br, const double *bi, + double *cr, double *ci + ) +{ + const double *pvalr, *pvali; + double *pcr=cr; + double *pci=ci; + int i,j; + int jj; + int rpntrb, rpntre; + int index, nvals; + + + for (i=0;i!=m;i++){ + *pcr++ = 0.0; + *pci++ = 0.0; + } + + pvalr = valr; + pvali = vali; + for (j=0;j!=k;j++){ + rpntrb = pntrb[j]; + rpntre = pntrb[j+1]; + for (jj=rpntrb;jj!=rpntre;jj++) { + index = indx[jj]; + if ( index == j ) { + cr[j] += (br[j] * (*pvalr)) - (bi[j] * (*pvali)); + ci[j] += (br[j] * (*pvali++)) + (bi[j] * (*pvalr++)); + continue; + } + if ( index > j ) { + cr[index] += (br[j] * (*pvalr)) - (bi[j] * (*pvali)); + ci[index] += (br[j] * (*pvali)) + (bi[j] * (*pvalr)); + + cr[j] += (br[index] * (*pvalr)) - (bi[index] * (*pvali)); + ci[j] += (br[index] * (*pvali++)) + (bi[index] * (*pvalr++)); + } + else { + pvalr++; + pvali++; + } + + } + } +} + + +/*C<-A\B; with Lower triangular A*/ +void CSC_VecTriangSlvLD_CAB_double( + const int m, + const double *val, + const int *indx, const int *pntrb, + const double *b, + double *c) +{ + int i, j, jb, je; + double *pc=c; + double z; + + for (i=0;i!=m;i++){ + *pc = b[i]; + pc++; + } + + pc=c; + for (i=0;i!=m;i++) { + jb = pntrb[i]; + je = pntrb[i+1]; + z = pc[i] / val[jb]; + pc[i] = z; + for (j=jb+1;j1) { + mexErrMsgTxt("Too many output arguments"); + } + + /* get edgel information */ + nr = mxGetM(in[0]); + nc = mxGetN(in[0]); + if ( nr*nc ==0 || nr != mxGetM(in[1]) || nc != mxGetN(in[1]) ) { + mexErrMsgTxt("Edge magnitude and phase shall be of the same image size"); + } + emag = mxGetPr(in[0]); + ephase = mxGetPr(in[1]); + np = nr * nc; + + /* get new index pair */ + if (!mxIsUint32(in[2]) | !mxIsUint32(in[3])) { + mexErrMsgTxt("Index pair shall be of type UINT32"); + } + if (mxGetM(in[3]) * mxGetN(in[3]) != np + 1) { + mexErrMsgTxt("Wrong index representation"); + } + pi = mxGetData(in[2]); + pj = mxGetData(in[3]); + + /* create output */ + out[0] = mxCreateSparse(np,np,pj[np],mxREAL); + if (out[0]==NULL) { + mexErrMsgTxt("Not enough memory for the output matrix"); + } + w = mxGetPr(out[0]); + ir = mxGetIr(out[0]); + jc = mxGetJc(out[0]); + + /* find my sigma */ + if (nargin<5) { + sigma = 0; + for (k=0; ksigma) { sigma = emag[k]; } + } + sigma = sigma / 6; + /* printf("sigma = %6.5f",sigma); */ + } else { + sigma = mxGetScalar(in[4]); + } + a = 0.5 / (sigma * sigma); + + /* computation */ + total = 0; + for (j=0; j= abs(dj)) { + slope = dj / di; + step = (iy>=jy) ? 1 : -1; + + iip1 = jy; + jjp1 = jx; + + + for (ii=0;ii maxori){ + maxori = z; + } + } + + iip1 = iip2; + jjp1 = jjp2; + phase1 = phase2; + } + + /* sample in j direction */ + } else { + slope = di / dj; + step = (ix>=jx) ? 1: -1; + + jjp1 = jx; + iip1 = jy; + + + for (jj=0;jj maxori){ + maxori = z; + } + + } + + iip1 = iip2; + jjp1 = jjp2; + phase1 = phase2; + } + } + + maxori = 0.5 * maxori; + maxori = exp(-maxori * maxori * a); + } + ir[total] = i; + + w[total] = maxori; + total = total + 1; + + } /* i */ + } /* j */ + + jc[np] = total; +} diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.dll b/SD-VBS/common/toolbox/MultiNcut/affinityic.dll new file mode 100755 index 0000000..67e4d64 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/affinityic.dll differ diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64 b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64 new file mode 100755 index 0000000..2648387 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexa64 differ diff --git a/SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx new file mode 100755 index 0000000..c296845 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/affinityic.mexglx differ diff --git a/SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m b/SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m new file mode 100755 index 0000000..d4dcb3c --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/batch_MNcut.m @@ -0,0 +1,48 @@ +data_dir = '/data/insecure/qihuizhu/baseball/Gray/train/'; +save_dir = '/home/songgang/project/MultiNcut/batch_result_MNcut'; + +num_segs = [20]; +imageSize = 200; + + +filelist = dir(fullfile(data_dir, '*.tif')); + +nb_file = max(size(filelist)); + + +tic; +for ii = 1:nb_file + fprintf(2, 'Segmenting image: %s ...\n', filelist(ii).name); + + img_filename = fullfile(data_dir, filelist(ii).name); + I=readimage(img_filename,imageSize); + + + [SegLabel,eigenVectors,eigenValues]= MNcut(I,num_segs); + + for j=1:size(SegLabel,3), + [gx,gy] = gradient(SegLabel(:,:,j)); + bw = (abs(gx)>0.1) + (abs(gy) > 0.1); + + figure(1);clf; J1=showmask(double(I),bw); imagesc(J1);axis image; axis off; + set(gca, 'Position', [0 0 1 1]); + set(gca, 'Position', [0 0 1 1]); + [PATHSTR,NAME,EXT,VERSN] = fileparts(filelist(ii).name); + print('-f1', '-djpeg90', fullfile(save_dir, sprintf('%s%s-%d.jpg', NAME,'-out', num_segs(j)))); + + + % cm = sprintf('print -djpeg %s/file%.4d-%.2d.jpg',OutputDir,id,num_segs(j)); disp(cm);eval(cm); + + + % figure(10);imagesc(SegLabel(:,:,j));axis image; axis off; + % set(gca, 'Position', [0 0 1 1]); + % cm = sprintf('print -djpeg %s/Seg%.4d-%.2d.jpg',OutputDir,id,num_segs(j));disp(cm);eval(cm); + +% keyboard; + end + + + +end; +toc; +fprintf(2, ' %d files done\n', nb_file); diff --git a/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c new file mode 100755 index 0000000..44af715 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/cimgnbmap.c @@ -0,0 +1,198 @@ +/*================================================================ +* function [i,j] = cimgnbmap([nr,nc], nb_r, sample_rate) +* computes the neighbourhood index matrix of an image, +* with each neighbourhood sampled. +* Input: +* [nr,nc] = image size +* nb_r = neighbourhood radius, could be [r_i,r_j] for i,j +* sample_rate = sampling rate, default = 1 +* Output: +* [i,j] = each is a column vector, give indices of neighbour pairs +* UINT32 type +* i is of total length of valid elements, 0 for first row +* j is of length nr * nc + 1 +* +* See also: imgnbmap.c, id2cind.m +* +* Examples: +* [i,j] = imgnbmap(10, 20); % [10,10] are assumed +* +* Stella X. Yu, Nov 12, 2001. + +% test sequence: +nr = 15; +nc = 15; +nbr = 1; +[i,j] = cimgnbmap([nr,nc], nbr); +mask = csparse(i,j,ones(length(i),1),nr*nc); +show_dist_w(rand(nr,nc),mask) + +*=================================================================*/ + +# include "mex.h" +# include "math.h" + +void mexFunction( + int nargout, + mxArray *out[], + int nargin, + const mxArray *in[] +) +{ + /* declare variables */ + int nr, nc, np, nb, total; + double *dim, sample_rate; + int r_i, r_j, a1, a2, b1, b2, self, neighbor; + int i, j, k, s, t, nsamp, th_rand, no_sample; + /* unsigned long *p, *qi, *qj; */ + unsigned int *p, *qi, *qj; + + /* check argument */ + if (nargin < 2) { + mexErrMsgTxt("Two input arguments required"); + } + if (nargout> 2) { + mexErrMsgTxt("Too many output arguments."); + } + + /* get image size */ + i = mxGetM(in[0]); + j = mxGetN(in[0]); + dim = mxGetData(in[0]); + nr = (int)dim[0]; + if (j>1 || i>1) { + nc = (int)dim[1]; + } else { + nc = nr; + } + np = nr * nc; + + /* get neighbourhood size */ + i = mxGetM(in[1]); + j = mxGetN(in[1]); + dim = mxGetData(in[1]); + r_i = (int)dim[0]; + if (j>1 || i>1) { + r_j = (int)dim[1]; + } else { + r_j = r_i; + } + if (r_i<0) { r_i = 0; } + if (r_j<0) { r_j = 0; } + + /* get sample rate */ + if (nargin==3) { + sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]); + } else { + sample_rate = 1; + } + /* prepare for random number generator */ + if (sample_rate<1) { + srand( (unsigned)time( NULL ) ); + th_rand = (int)ceil((double)RAND_MAX * sample_rate); + no_sample = 0; + } else { + sample_rate = 1; + th_rand = RAND_MAX; + no_sample = 1; + } + + /* figure out neighbourhood size */ + + nb = (r_i + r_i + 1) * (r_j + r_j + 1); + if (nb>np) { + nb = np; + } + nb = (int)ceil((double)nb * sample_rate); +/*printf("nb=%d\n",nb);*/ + /* intermediate data structure */ + /* p = mxCalloc(np * (nb+1), sizeof(unsigned long)); */ + p = mxCalloc(np * (nb+1), sizeof(unsigned int)); + if (p==NULL) { + mexErrMsgTxt("Not enough space for my computation."); + } + + /* computation */ + total = 0; + for (j=0; j=nc) { b2 = nc-1; } + + /* i range */ + a1 = i - r_i; + if (a1<0) { a1 = 0; } + a2 = i + r_i; + if (a2>=nr) { a2 = nr-1; } + + /* number of more samples needed */ + nsamp = nb - p[self]; + /*if (nsamp<0) + {printf("nsamp=%d\n",nsamp);}*/ + k = 0; + t = b1; + s = i + 1; + if (s>a2) { + s = a1; + t = t + 1; + } + + + while (ka2) { + s = a1; + t = t + 1; + } + } /* k */ + total = total + p[self]; + } /* i */ + + } /* j */ + + /* i, j */ + + out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL); + out[1] = mxCreateNumericMatrix(np+1, 1, mxUINT32_CLASS, mxREAL); + qi = mxGetData(out[0]); + qj = mxGetData(out[1]); + + if (out[0]==NULL || out[1]==NULL) { + mexErrMsgTxt("Not enough space for the output matrix."); + } + + total = 0; + for (j=0; j 2) { + mexErrMsgTxt("Too many output arguments."); + } + + /* get image size */ + i = mxGetM(in[0]); + j = mxGetN(in[0]); + dim = mxGetData(in[0]); + nr = (int)dim[0]; + if (j>1 || i>1) { + nc = (int)dim[1]; + } else { + nc = nr; + } + np = nr * nc; + + /* get neighbourhood size */ + i = mxGetM(in[1]); + j = mxGetN(in[1]); + dim = mxGetData(in[1]); + r_i = (int)dim[0]; + if (j>1 || i>1) { + r_j = (int)dim[1]; + } else { + r_j = r_i; + } + if (r_i<0) { r_i = 0; } + if (r_j<0) { r_j = 0; } + + /* get sample rate */ + if (nargin==3) { + sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]); + } else { + sample_rate = 1; + } + /* prepare for random number generator */ + if (sample_rate<1) { + srand( (unsigned)time( NULL ) ); + th_rand = (int)ceil((double)RAND_MAX * sample_rate); + no_sample = 0; + } else { + sample_rate = 1; + th_rand = RAND_MAX; + no_sample = 1; + } + + /* figure out neighbourhood size */ + + nb = (r_i + r_i) * (r_j + r_j)+1; + if (nb>np) { + nb = np; + } + nb = (int)ceil((double)nb * sample_rate); +/*printf("nb=%d\n",nb);*/ + /* intermediate data structure */ + /* p = mxCalloc(np * (nb+1), sizeof(unsigned long)); */ + p = mxCalloc(np * (nb+1), sizeof(unsigned int)); + if (p==NULL) { + mexErrMsgTxt("Not enough space for my computation."); + } + + /* computation */ + total = 0; + for (j=0; j=nc) { b2 = nc-1; } + + /* i range */ + /*a1 = i - r_i; + if (a1<0) { a1 = 0; }*/ + a2 = i + r_i; + if (a2>=nr) { a2 = nr-1; } + + /* number of more samples needed */ + nsamp = nb - p[self]; + /*if (nsamp<0) + {printf("nsamp=%d\n",nsamp);}*/ + k = 0; + t = b1; + s = i + 1; + if (s>a2) { + s = i; + t = t + 1; + } + + + while (ka2) { s = i; t = t + 1; + } + } + else {t=t+1;} + } /* k */ + total = total + p[self]; + } /* i */ + + } /* j */ + + /* i, j */ + + out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL); + out[1] = mxCreateNumericMatrix(np+1, 1, mxUINT32_CLASS, mxREAL); + qi = mxGetData(out[0]); + qj = mxGetData(out[1]); + + if (out[0]==NULL || out[1]==NULL) { + mexErrMsgTxt("Not enough space for the output matrix."); + } + + total = 0; + for (j=0; j 2) { + mexErrMsgTxt("Too many output arguments."); + } + + + /* get image size */ + + i = mxGetM(in[0]); + j = mxGetN(in[0]); + dim = mxGetData(in[0]); + nr = (int)dim[0]; + if (j>1 || i>1) { + nc = (int)dim[1]; + } else { + nc = nr; + } + np = nr * nc; + + + /* get neighbourhood size */ + i = mxGetM(in[1]); + j = mxGetN(in[1]); + dim = mxGetData(in[1]); + r_i = (int)dim[0]; + + if (j>1 || i>1) { + r_j = (int)dim[1]; + } else { + r_j = r_i; + } + + if (r_i<0) { r_i = 0; } + + if (r_j<0) { r_j = 0; } + + + + /* get sample rate */ + + if (nargin==3) { + + sample_rate = (mxGetM(in[2])==0) ? 1: mxGetScalar(in[2]); + + } else { + + sample_rate = 1; + + } + + /* prepare for random number generator */ + if (sample_rate<1) { + srand( (unsigned)time( NULL ) ); + + th_rand = (int)ceil((double)RAND_MAX * sample_rate); + no_sample = 0; + } else { + + sample_rate = 1; + th_rand = RAND_MAX; + no_sample = 1; + } + + + /* figure out neighbourhood size */ + + nb = (4*r_i) + (4*r_j)+1; + if (nb>np) { + nb = np; + } + nb = (int)ceil((double)nb * sample_rate); + +/*printf("nb=%d\n",nb);*/ + + /* intermediate data structure */ + + /* p = mxCalloc(np * (nb+1), sizeof(unsigned long));*/ + p = mxCalloc(np * (nb+1), sizeof(unsigned int)); + + if (p==NULL) { + + mexErrMsgTxt("Not enough space for my computation."); + + } + + + + /* computation */ + + total = 0; + for (j=0; j=nc) { b2 = nc-1; } + + + /* i range */ + /*a1 = i - r_i; + + if (a1<0) { a1 = 0; }*/ + a2 = i + r_i; + if (a2>=nr) { a2 = nr-1; } + + + /* number of more samples needed */ + + nsamp = nb - p[self]; + + /*if (nsamp<0) + {printf("nsamp=%d\n",nsamp);}*/ + k = 0; + t = b1; + s = i + 1; + + if (s>a2) { + s = i; + t = t + 1; + } + + + + while (ka2) { + t = t + 1; + if (i+j-t>=0) + {s = i+j-t;} + else + {s=i;} + } + } + else { + if (s==i+j-t) {s=i;} + else{ if (s==i && s+t-j=0) + {s = i+j-t;} + else + {s=i;} + } + + } + } + + } /* k */ + + total = total + p[self]; + } /* i */ + + } /* j */ + + + + /* i, j */ + + out[0] = mxCreateNumericMatrix(total, 1, mxUINT32_CLASS, mxREAL); + + out[1] = mxCreateNumericMatrix(np+1, 1, mxUINT32_CLASS, mxREAL); + + qi = mxGetData(out[0]); + + qj = mxGetData(out[1]); + + + if (out[0]==NULL || out[1]==NULL) { + + mexErrMsgTxt("Not enough space for the output matrix."); + + } + + + + total = 0; + + for (j=0; j1) + for j=i-1:-1:1 + Wpp{i}=[C{j,i}',Wpp{i}]; + end + end + if (i1) +% for j=i-1:-1:1 +% Wpp{i}=[sparse(p(i)*q(i),p(j)*q(j)),Wpp{i}]; +% end +% end +% if (i2 + for i=3:n + piqi=p(i)*q(i); + if i~=n + constraintMat=[constraintMat,[sparse(sum(p(1:i-2).*q(1:i-2)),piqi);-C{i-1,i};speye(piqi);sparse(pq-sum(p(2:i).*q(2:i)),piqi)]]; + else + constraintMat=[constraintMat,[sparse(sum(p(1:i-2).*q(1:i-2)),piqi);-C{i-1,i};speye(piqi)]]; + end + end + end + + % saving useful information + %subgrids, p and q + data.subgrid=subgrid; + data.p=p; + data.q=q; diff --git a/SD-VBS/common/toolbox/MultiNcut/discretisation.m b/SD-VBS/common/toolbox/MultiNcut/discretisation.m new file mode 100755 index 0000000..70b5650 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/discretisation.m @@ -0,0 +1,49 @@ +function [SegLabel,EigenVectors]=discretisation(EigenVectors,nr,nc) +% +% EigenvectorsDiscrete=discretisation(EigenVectors) +% +% Input: EigenVectors = continuous Ncut vector, size = ndata x nbEigenvectors +% Output EigenvectorsDiscrete = discrete Ncut vector, size = ndata x nbEigenvectors +% +% Timothee Cour, Stella Yu, Jianbo Shi, 2004 + +[n,k]=size(EigenVectors); + +vm = sqrt(sum(EigenVectors.*EigenVectors,2)); +EigenVectors = EigenVectors./repmat(vm,1,k); + +R=zeros(k); +R(:,1)=EigenVectors(1+round(rand(1)*(n-1)),:)'; +c=zeros(n,1); +for j=2:k + c=c+abs(EigenVectors*R(:,j-1)); + [minimum,i]=min(c); + R(:,j)=EigenVectors(i,:)'; +end + +lastObjectiveValue=0; +exitLoop=0; +nbIterationsDiscretisation = 0; +nbIterationsDiscretisationMax = 20;%voir +while exitLoop== 0 + nbIterationsDiscretisation = nbIterationsDiscretisation + 1 ; + EigenvectorsDiscrete = discretisationEigenVectorData(EigenVectors*R); + [U,S,V] = svd(EigenvectorsDiscrete'*EigenVectors,0); + NcutValue=2*(n-trace(S)); + + if abs(NcutValue-lastObjectiveValue) < eps | nbIterationsDiscretisation > nbIterationsDiscretisationMax + exitLoop=1; + else + lastObjectiveValue = NcutValue; + R=V*U'; + end +end + +%%%% + +SegLabel = zeros(nr,nc); +for j=1:size(EigenvectorsDiscrete,2), + SegLabel = SegLabel + j*reshape(EigenvectorsDiscrete(:,j),nr,nc); +end +EigenVectors = reshape(EigenVectors,nr,nc,size(EigenVectors,2)); + diff --git a/SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m b/SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m new file mode 100755 index 0000000..4626e3d --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/discretisationEigenVectorData.m @@ -0,0 +1,12 @@ +function Y = discretisationEigenVectorData(EigenVector) +% Y = discretisationEigenVectorData(EigenVector) +% +% discretizes previously rotated eigenvectors in discretisation +% Timothee Cour, Stella Yu, Jianbo Shi, 2004 + +[n,k]=size(EigenVector); + + +[Maximum,J]=max(EigenVector'); + +Y=sparse(1:n,J',1,n,k); diff --git a/SD-VBS/common/toolbox/MultiNcut/doog1.m b/SD-VBS/common/toolbox/MultiNcut/doog1.m new file mode 100755 index 0000000..dd8e87b --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/doog1.m @@ -0,0 +1,32 @@ +function H=doog1(sig,r,th,N); +% H=doog1(sig,r,th,N); + + +% by Serge Belongie + +no_pts=N; % no. of points in x,y grid + +[x,y]=meshgrid(-(N/2)+1/2:(N/2)-1/2,-(N/2)+1/2:(N/2)-1/2); + +phi=pi*th/180; +sigy=sig; +sigx=r*sig; +R=[cos(phi) -sin(phi); sin(phi) cos(phi)]; +C=R*diag([sigx,sigy])*R'; + +X=[x(:) y(:)]; + +Gb=gaussian(X,[0 0]',C); +Gb=reshape(Gb,N,N); + +m=R*[0 sig]'; + +a=1; +b=-1; + +% make odd-symmetric filter +Ga=gaussian(X,m/2,C); +Ga=reshape(Ga,N,N); +Gb=rot90(Ga,2); +H=a*Ga+b*Gb; + diff --git a/SD-VBS/common/toolbox/MultiNcut/doog2.m b/SD-VBS/common/toolbox/MultiNcut/doog2.m new file mode 100755 index 0000000..a0511cb --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/doog2.m @@ -0,0 +1,38 @@ +function G=doog2(sig,r,th,N); +% G=doog2(sig,r,th,N); +% Make difference of offset gaussians kernel +% theta is in degrees +% (see Malik & Perona, J. Opt. Soc. Amer., 1990) +% +% Example: +% >> imagesc(doog2(1,12,0,64,1)) +% >> colormap(gray) + +% by Serge Belongie + +no_pts=N; % no. of points in x,y grid + +[x,y]=meshgrid(-(N/2)+1/2:(N/2)-1/2,-(N/2)+1/2:(N/2)-1/2); + +phi=pi*th/180; +sigy=sig; +sigx=r*sig; +R=[cos(phi) -sin(phi); sin(phi) cos(phi)]; +C=R*diag([sigx,sigy])*R'; + +X=[x(:) y(:)]; + +Gb=gaussian(X,[0 0]',C); +Gb=reshape(Gb,N,N); + +m=R*[0 sig]'; +Ga=gaussian(X,m,C); +Ga=reshape(Ga,N,N); +Gc=rot90(Ga,2); + +a=-1; +b=2; +c=-1; + +G = a*Ga + b*Gb + c*Gc; + diff --git a/SD-VBS/common/toolbox/MultiNcut/eigSolve.m b/SD-VBS/common/toolbox/MultiNcut/eigSolve.m new file mode 100755 index 0000000..cfb64fc --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/eigSolve.m @@ -0,0 +1,5 @@ +function [eigenVectors,s]= eigSolve (multiWpp,ConstraintMat,data) + +[v,s] = quickNcutHardBiais2(multiWpp,ConstraintMat,data.dataGraphCut.nbEigenValues,data.dataGraphCut); +v=v(1:data.p(1)*data.q(1),:); +eigenVectors=reshape(v,data.p(1),data.q(1),data.dataGraphCut.nbEigenValues); diff --git a/SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m b/SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m new file mode 100755 index 0000000..9c84e96 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/fft_filt_2.m @@ -0,0 +1,29 @@ +function FI=fft_filt_2(V,FB,sf); +% FI=fft_filt_2(V,FB,sf); +% fft-based filtering +% requires image to be called "V" +% and filters to be in FB +% sf is the subsampling factor +% +% FI is the result + +[M1,M2,N3]=size(FB); +% prepare FFT of image for filtering +[N1,N2]=size(V); +I=zeros(size(V)+[M1-1 M2-1]); +I(1:N1,1:N2)=V; +N1s=length(1:sf:N1); +N2s=length(1:sf:N2); +IF=fft2(I); +FI=zeros(N1s,N2s,N3); + +% apply filters +for n=1:N3; + f=rot90(FB(:,:,n),2); + fF=fft2(f,N1+M1-1,N2+M2-1); + IfF=IF.*fF; + If=real(ifft2(IfF)); + If=If(ceil((M1+1)/2):ceil((M1+1)/2)+N1-1,ceil((M2+1)/2):ceil((M2+1)/2)+N2-1); + FI(:,:,n)=If(1:sf:N1,1:sf:N2); +end + diff --git a/SD-VBS/common/toolbox/MultiNcut/gaussian.m b/SD-VBS/common/toolbox/MultiNcut/gaussian.m new file mode 100755 index 0000000..509b129 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/gaussian.m @@ -0,0 +1,31 @@ +function p=gaussian(x,m,C); +% p=gaussian(x,m,C); +% +% Evaluate the multi-variate density with mean vector m and covariance +% matrix C for the input vector x. +% +% p=gaussian(X,m,C); +% +% Vectorized version: Here X is a matrix of column vectors, and p is +% a vector of probabilities for each vector. + +d=length(m); + +if size(x,1)~=d + x=x'; +end +N=size(x,2); + +detC = det(C); +if rcond(C)0) { + ir[k] = mxGetIr(in[k]); + jc[k] = mxGetJc(in[k]); + } + } + + out[0] = mxCreateDoubleMatrix(m[1],1,mxREAL); + y = mxGetPr(out[0]); + + + y1 = mxCalloc(n[2], sizeof(double)); + y2 = mxCalloc(m[3], sizeof(double)); + y2bis = mxCalloc(m[3], sizeof(double)); + y3 = mxCalloc(m[1], sizeof(double)); + y4 = mxCalloc(m[1], sizeof(double)); + y5 = mxCalloc(n[2], sizeof(double)); + y5bis = mxCalloc(n[2], sizeof(double)); + y6 = mxCalloc(n[2], sizeof(double)); + + CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],pr[0],y1); + CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y1,y2bis); + CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y2bis,y2); + for(k=0;k0) { + ir[k] = mxGetIr(in[k]); + jc[k] = mxGetJc(in[k]); + } + } + + out[0] = mxCreateDoubleMatrix(m[1],1,mxREAL); + y = mxGetPr(out[0]); + + + y1 = mxCalloc(n[2], sizeof(double)); + y2 = mxCalloc(m[3], sizeof(double)); + y2bis = mxCalloc(m[3], sizeof(double)); + y3 = mxCalloc(m[1], sizeof(double)); + y4 = mxCalloc(m[1], sizeof(double)); + y5 = mxCalloc(n[2], sizeof(double)); + y5bis = mxCalloc(n[2], sizeof(double)); + y6 = mxCalloc(n[2], sizeof(double)); + + CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],pr[0],y1); + CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y1,y2bis); + CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y2bis,y2); + for(k=0;k1) { + mexErrMsgTxt("Too many output arguments"); + } + + /* get edgel information */ + nr = mxGetM(in[0]); + nc = mxGetN(in[0]); + if ( nr*nc ==0 || nr != mxGetM(in[1]) || nc != mxGetN(in[1]) ) { + mexErrMsgTxt("Edge magnitude and phase shall be of the same image size"); + } + emag = mxGetPr(in[0]); + ephase = mxGetPr(in[1]); + np = nr * nc; + + /*get subgrid information*/ + + tmp = mxGetData(in[5]); + nrSubgrid = (int)tmp[0]; + tmp = mxGetData(in[6]); + ncSubgrid = (int)tmp[0]; + + /* printf("nrSubgrid=%d\n",nrSubgrid); + printf("ncSubgrid=%d\n",ncSubgrid); */ + if (nrSubgrid* ncSubgrid != mxGetM(in[4])*mxGetN(in[4])) { + mexErrMsgTxt("Error in the size of the subgrid"); + } + subgrid = mxGetData(in[4]); + nW = nrSubgrid * ncSubgrid; + + /* get new index pair */ + if (!mxIsUint32(in[2]) | !mxIsUint32(in[3])) { + mexErrMsgTxt("Index pair shall be of type UINT32"); + } + if (mxGetM(in[3]) * mxGetN(in[3]) != nW + 1) { + mexErrMsgTxt("Wrong index representation"); + } + pi = mxGetData(in[2]); + pj = mxGetData(in[3]); + subpi = mxGetData(in[7]); + /*{printf("pi[50] = %d\n",pi[50]);} + {printf("subpi[5] = %d\n",subpi[5]);} + {printf("subpi[6] = %d\n",subpi[6]);} + {printf("subpi[4] = %d\n",subpi[4]);}*/ + + /* create output */ /*!!!!!!!!!!!!!!!!!!!!!!!changer taille output!!!!!!!!!!*/ + out[0] = mxCreateSparse(nW,nW,pj[nW],mxREAL); + if (out[0]==NULL) { + mexErrMsgTxt("Not enough memory for the output matrix"); + } + w = mxGetPr(out[0]); + ir = mxGetIr(out[0]); + jc = mxGetJc(out[0]); + + /* find my sigma */ + if (nargin<9) { + sigma = 0; + for (k=0; ksigma) { sigma = emag[k]; } + } + sigma = sigma / 6; + printf("sigma = %6.5f",sigma); + } else { + sigma = mxGetScalar(in[8]); + } + a = 0.5 / (sigma * sigma); + + /* computation */ + total = 0; + + for (j=0; j= abs(dj)) { + slope = dj / di; + step = (iy>=jy) ? 1 : -1; + + iip1 = jy; + jjp1 = jx; + + + for (ii=0;ii maxori){ + maxori = z; + } + } + + iip1 = iip2; + jjp1 = jjp2; + phase1 = phase2; + } + + /* sample in j direction */ + } else { + slope = di / dj; + step = (ix>=jx) ? 1: -1; + + jjp1 = jx; + iip1 = jy; + + + for (jj=0;jj maxori){ + maxori = z; + } + + } + + iip1 = iip2; + jjp1 = jjp2; + phase1 = phase2; + } + } + + maxori = 0.5 * maxori; + maxori = exp(-maxori * maxori * a); + } + /*if (total<20) {printf("subpi[k] = %d\n",subpi[k]);}*/ + + ir[total] = (int)subpi[k]; +/*if (total<20) {printf("ir[total] = %d\n",ir[total]);}*/ + w[total] = maxori; + total = total + 1; + + } /* i */ + } /*j*/ + /*printf("total = %d\n",total);*/ + /*printf("ir[100] = %d\n",ir[100]);*/ + jc[nW] = total; +} diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll new file mode 100755 index 0000000..f9b6bc9 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.dll differ diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64 b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64 new file mode 100755 index 0000000..9a2a3cc Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexa64 differ diff --git a/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx new file mode 100755 index 0000000..ec11a0e Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/multiAffinityic.mexglx differ diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c new file mode 100755 index 0000000..0cae523 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/multiIntensityFirstLayer.c @@ -0,0 +1,126 @@ +/*================================================================ +* function w = intensityWppc(image,pi,pj,rMin,sigmaX,sigmaIntensite,valeurMinW ) +* Input: +* [pi,pj] = index pair representation for MALTAB sparse matrices +* Output: +* w = affinity with IC at [pi,pj] +* + +pixels i and j (corresponding to the sampling in pi,pj) are fully connected when d(i,j) <= rmin; + +% test sequence +f = synimg(10); +[i,j] = cimgnbmap(size(f),2); +[ex,ey,egx,egy] = quadedgep(f); +a = affinityic(ex,ey,egx,egy,i,j) +show_dist_w(f,a); + +* Stella X. Yu, Nov 19, 2001. +*=================================================================*/ + +# include "mex.h" +# include "math.h" + +void mexFunction( + int nargout, + mxArray *out[], + int nargin, + const mxArray *in[] +) +{ + /* declare variables */ + int nr, nc, np, total; + int i, j, k, ix, iy, jx, jy; + int *ir, *jc; + int squareDistance; + /* unsigned long *pi, *pj; */ +unsigned int *pi, *pj; + double *w; + + double temp,a1,a2,wij; + double rMin; + double sigmaX, sigmaIntensite,valeurMinW; + double *image; + + /* check argument */ + if (nargin<7) { + mexErrMsgTxt("Four input arguments required"); + } + if (nargout>1) { + mexErrMsgTxt("Too many output arguments"); + } + + /* get edgel information */ + nr = mxGetM(in[0]); + nc = mxGetN(in[0]); + np = nr * nc; + + image = mxGetPr(in[0]); + + + + /* get new index pair */ + if (!mxIsUint32(in[1]) | !mxIsUint32(in[2])) { + mexErrMsgTxt("Index pair shall be of type UINT32"); + } + if (mxGetM(in[2]) * mxGetN(in[2]) != np + 1) { + mexErrMsgTxt("Wrong index representation"); + } + pi = mxGetData(in[1]); + pj = mxGetData(in[2]); + + /* create output */ + out[0] = mxCreateSparse(np,np,pj[np],mxREAL); + if (out[0]==NULL) { + mexErrMsgTxt("Not enough memory for the output matrix"); + } + w = mxGetPr(out[0]); + ir = mxGetIr(out[0]); + jc = mxGetJc(out[0]); + + rMin = mxGetScalar(in[3]); + sigmaX = mxGetScalar(in[4]); + sigmaIntensite= mxGetScalar(in[5]); + valeurMinW = mxGetScalar(in[6]); + + a1 = 1.0/ (sigmaX*sigmaX); + a2 = 1.0 / (sigmaIntensite*sigmaIntensite ); + + /* computation */ + total = 0; + for (j=0; j1) { + mexErrMsgTxt("Too many output arguments"); + } + + /* get edgel information */ + nr = mxGetM(in[0]); + nc = mxGetN(in[0]); + np = nr * nc; + /*printf("size: %d, %d, %d\n", nc, nr, np); */ + image = mxGetPr(in[0]); + + + /*get subgrid information*/ + + tmp = mxGetData(in[8]); + nrSubgrid = (int)tmp[0]; + + /* printf("image end = %f ", image[np-1]); */ + + tmp = mxGetData(in[9]); + ncSubgrid = (int)tmp[0]; + + if (nrSubgrid* ncSubgrid != mxGetM(in[7])*mxGetN(in[7])) { + mexErrMsgTxt("Error in the size of the subgrid"); + } + subgrid = mxGetData(in[7]); + nW = nrSubgrid * ncSubgrid; + + + + + /* get new index pair */ + if (!mxIsUint32(in[1]) | !mxIsUint32(in[2])) { + mexErrMsgTxt("Index pair shall be of type UINT32"); + } + if (mxGetM(in[2]) * mxGetN(in[2]) != nW + 1) { + mexErrMsgTxt("Wrong index representation"); + } + pi = mxGetData(in[1]); + pj = mxGetData(in[2]); + subpi = mxGetData(in[10]); + + /* create output */ + out[0] = mxCreateSparse(nW,nW,pj[nW],mxREAL); + if (out[0]==NULL) { + mexErrMsgTxt("Not enough memory for the output matrix"); + } + + w = mxGetPr(out[0]); + ir = mxGetIr(out[0]); + jc = mxGetJc(out[0]); + + + rMin = mxGetScalar(in[3]); + sigmaX = mxGetScalar(in[4]); + sigmaIntensite= mxGetScalar(in[5]); + valeurMinW = mxGetScalar(in[6]); + + a1 = 1.0/ (sigmaX*sigmaX); + a2 = 1.0 / (sigmaIntensite*sigmaIntensite ); + + + + /* computation */ + total = 0; + for (j=0; jnp-1)) {printf("badddddd!");} + /* printf("t = %d\n",t); + printf("j = %d\n",j); */ + jc[j] = total; + jx = t / nr; /* col */ + jy = t % nr; /* row */ + /* printf("pj[j+1] = %d\n",pj[j+1]); */ + + for (k=pj[j]; k5000*5000 ) {printf("trouble! [%d,%d]\n",k,(int)subpi[k]);} */ + + w[total] = wij; + + total = total + 1; + + } /* i */ + + + } /* j */ + + jc[nW] = total; +} diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll new file mode 100755 index 0000000..16146c1 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.dll differ diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64 b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64 new file mode 100755 index 0000000..525ca62 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexa64 differ diff --git a/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx new file mode 100755 index 0000000..a17fa03 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/multiIntensityWppc.mexglx differ diff --git a/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m b/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m new file mode 100755 index 0000000..5041377 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m @@ -0,0 +1,188 @@ +% function [xs,ys,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold); +% Input: +% I = image +% par = vector for 4 parameters +% [number of filter orientations, number of scale, filter size, elongation] +% To use default values, put 0. +% threshold = threshold on edge strength +% Output: +% [x,y,gx,gy] = locations and gradients of an ordered list of edgels +% x,y could be horizontal or vertical or 45 between pixel sites +% but it is guaranteed that there [floor(y) + (floor(x)-1)*nr] +% is ordered and unique. In other words, each edgel has a unique pixel id. +% par = actual par used +% threshold = actual threshold used +% mag = edge magnitude +% mage = phase map +% g = gradient map at each pixel +% [FIe,FIo] = odd and even filter outputs +% mago = odd filter output of optimum orientation + +% Stella X. Yu, 2001 + +% This is the multi scale version of the filtering +% For the moment the parameters are defined by default. See line 34 + + +function [x,y,gx,gy,par,threshold,mag_s,mage,g,FIe,FIo,mago] = quadedgep2(I,par,data,threshold); + + +if nargin<4 | isempty(threshold), + threshold = 0.1; +end + +[r,c] = size(I); +def_par = [4,30,3]; + +display_on = 1; + +% take care of parameters, any missing value is substituted by a default value +if nargin<2 | isempty(par), + par = def_par; +end +% par(end+1:4)=0; +% par = par(:); +% j = (par>0); +% have_value = [ j, 1-j ]; +% j = 1; n_filter = have_value(j,:) * [par(j); def_par(j)]; +% j = 2; n_scale = have_value(j,:) * [par(j); def_par(j)]; +% j = 3; winsz = have_value(j,:) * [par(j); def_par(j)]; +% j = 4; enlong = have_value(j,:) * [par(j); def_par(j)]; + +n_ori = par(1); %if it doesn't work, par<-def_par + +winsz = par(2); +enlong = par(3); + +% always make filter size an odd number so that the results will not be skewed +j = winsz/2; +if not(j > fix(j) + 0.1), + winsz = winsz + 1; +end + +% filter the image with quadrature filters +if (isempty(data.W.scales)) + error ('no scales entered'); +end + +n_scale=length(data.W.scales); +filter_scales=data.W.scales; +% +% if strcmp(data.dataWpp.mode,'multiscale') +% n_scale=size((data.dataWpp.scales),2); +% filter_scales=data.dataWpp.scales; +% else +% filter_scales=data.dataWpp.scales(1); +% n_scale=1; +% end +% if n_scale>0&&strcmp(data.dataWpp.mode,'multiscale') +% if (~isempty(data.dataWpp.scales)) +% filter_scales=data.dataWpp.scales; +% else +% filter_scales=zeros(1,n_scale); +% +% for i=1:n_scale, +% filter_scales(i)=(sqrt(2))^(i-1); +% end +% data.dataWpp.scales=filter_scales; +% end +% else filter_scale=1; +% data.dataWpp.scales=filter_scales; +% end +% +% %%%%%%% juste pour que ca tourne +% if ~strcmp(data.dataWpp.mode,'multiscale') +% filter_scales=data.dataWpp.scales(1); +% n_scale=4; +% end +% %%%%%%%%%%%% + +FBo = make_filterbank_odd2(n_ori,filter_scales,winsz,enlong); +FBe = make_filterbank_even2(n_ori,filter_scales,winsz,enlong); +n = ceil(winsz/2); +f = [fliplr(I(:,2:n+1)), I, fliplr(I(:,c-n:c-1))]; +f = [flipud(f(2:n+1,:)); f; flipud(f(r-n:r-1,:))]; +FIo = fft_filt_2(f,FBo,1); +FIo = FIo(n+[1:r],n+[1:c],:); +FIe = fft_filt_2(f,FBe,1); +FIe = FIe(n+[1:r],n+[1:c],:); + +% compute the orientation energy and recover a smooth edge map +% pick up the maximum energy across scale and orientation +% even filter's output: as it is the second derivative, zero cross localize the edge +% odd filter's output: orientation + +[nr,nc,nb] = size(FIe); + +FIe = reshape(FIe, nr,nc,n_ori,length(filter_scales)); +FIo = reshape(FIo, nr,nc,n_ori,length(filter_scales)); + +mag_s = zeros(nr,nc,n_scale); +mag_a = zeros(nr,nc,n_ori); + +mage = zeros(nr,nc,n_scale); +mago = zeros(nr,nc,n_scale); +mage = zeros(nr,nc,n_scale); +mago = zeros(nr,nc,n_scale); + + + +for i = 1:n_scale, + mag_s(:,:,i) = sqrt(sum(FIo(:,:,:,i).^2,3)+sum(FIe(:,:,:,i).^2,3)); + mag_a = sqrt(FIo(:,:,:,i).^2+FIe(:,:,:,i).^2); + [tmp,max_id] = max(mag_a,[],3); + + base_size = nr * nc; + id = [1:base_size]'; + mage(:,:,i) = reshape(FIe(id+(max_id(:)-1)*base_size+(i-1)*base_size*n_ori),[nr,nc]); + mago(:,:,i) = reshape(FIo(id+(max_id(:)-1)*base_size+(i-1)*base_size*n_ori),[nr,nc]); + + mage(:,:,i) = (mage(:,:,i)>0) - (mage(:,:,i)<0); + + if display_on, + ori_incr=pi/n_ori; % to convert jshi's coords to conventional image xy + ori_offset=ori_incr/2; + theta = ori_offset+([1:n_ori]-1)*ori_incr; % orientation detectors + % [gx,gy] are image gradient in image xy coords, winner take all + + ori = theta(max_id); + ori = ori .* (mago(:,:,i)>0) + (ori + pi).*(mago(:,:,i)<0); + gy{i} = mag_s(:,:,i) .* cos(ori); + gx{i} = -mag_s(:,:,i) .* sin(ori); + g = cat(3,gx{i},gy{i}); + + % phase map: edges are where the phase changes + mag_th = max(max(mag_s(:,:,i))) * threshold; + eg = (mag_s(:,:,i)>mag_th); + h = eg & [(mage(2:r,:,i) ~= mage(1:r-1,:,i)); zeros(1,nc)]; + v = eg & [(mage(:,2:c,i) ~= mage(:,1:c-1,i)), zeros(nr,1)]; + [y{i},x{i}] = find(h | v); + k = y{i} + (x{i}-1) * nr; + h = h(k); + v = v(k); + y{i} = y{i} + h * 0.5; % i + x{i} = x{i} + v * 0.5; % j + t = h + v * nr; + gx{i} = g(k) + g(k+t); + k = k + (nr * nc); + gy{i} = g(k) + g(k+t); + +% figure(1); +% clf; +% imagesc(I);colormap(gray); +% hold on; +% quiver(x,y,gx,gy); hold off; +% title(sprintf('scale = %d, press return',i)); + + % pause; + 0; +else + x = []; + y = []; + gx = []; + gy =[]; + g= []; + end +end + + diff --git a/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m b/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m new file mode 100755 index 0000000..3ca1046 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/quickNcutHardBiais2.m @@ -0,0 +1,187 @@ +% function [v,s] = quickNcutHardBiais(W,U,nbEigenValues,dataGraphCut) +%ligne 35 : changement tim +%[v,s] = ncut(W,nbEigenValues,[],offset); +%devient : +%[v,s] = tim_ncut_rapide(W,nbEigenValues,[],offset); +%et eigs devient tim_eigs + +% Input: +% W = affinity matrix +% U = hard constraint matrix, could be a cell of partial grouping +% nbEigenValues = number of eigenvectors +% offset = regularization factor, default = 0 +% Output: +% v = eigenvector +% s = eigenvalue of (W,d), s.t. U' * x = 0. + +% call eigs using my own * operation + +% Stella X. Yu, Jan 2002. + +function [v,s] = quickNcutHardBiais2(W,U,nbEigenValues,dataGraphCut)%voir : rajouter sigma +n = size(W,1); +nbEigenValues = min(nbEigenValues,n); + +offset = dataGraphCut.offset; +%offset = 2; + +% degrees and regularization +d = sum(abs(W),2); +dr = 0.5 * (d - sum(W,2)); +d = d + offset * 2; +dr = dr + offset; +W = W + spdiags(dr,0,n,n); +clear dr + +% normalize +Dinvsqrt = 1./sqrt(d+eps); +P = spmtimesd(W,Dinvsqrt,Dinvsqrt); +clear W; + +% if max(max(P-P')) < 1e-10 %ou eps +% %S = sparse(1:n,1:n,0.5); +% P =max(P,P'); +% % P=S*(P+P'); +% %P=0.5*(P+P'); +% options.issym = 1; +% end +P = sparsifyc(P,dataGraphCut.valeurMin); +options.issym = 1; + +Ubar = spmtimesd(U,Dinvsqrt,[]); +%Ubar = sparsifyc(Ubar,dataGraphCut.valeurMin); %voir + +options.disp = dataGraphCut.verbose; +options.maxit = dataGraphCut.maxiterations; +options.tol = dataGraphCut.eigsErrorTolerance; + +options.v0 = ones(size(P,1),1);%voir + +options.p = max(35,2*nbEigenValues); %voir +options.p = min(options.p , n); + +% nouvelle idee : factorisation de Cholesky +C=Ubar'*Ubar; +%permutation = symamd(C); +%R = cholinc(C(permutation,permutation)); +t_chol_Ubar = cputime; +[R,ooo] = cholinc(C,'0'); +t_chol_Ubar = cputime - t_chol_Ubar; +%if error occurs, check if Ubar = sparsifyc(Ubar,dataGraphCut.valeurMin); +%sparsifies too much + + +% compute H = (Ubar'*Ubar)^(-1) +% t_inv_H = cputime; +% H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); %changer +% t_inv_H = cputime - t_inv_H; +% H = sparsifyc(H,dataGraphCut.valeurMin); +% tEigs = cputime; +% if options.issym & max(max(H-H')) < 1e-10 +% [vbar,s,convergence] = tim_eigs(@mex_projection_inv_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,triu(H)); +% else +% [vbar,s,convergence] = tim_eigs(@mex_projection_inv,n,nbEigenValues,'lm',options,P,Ubar,H); +% end +% tEigs = cputime - tEigs; +% + + + +R = sparsifyc(R,dataGraphCut.valeurMin); +tEigs = cputime; +if options.issym + [vbar,s,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,tril(P),Ubar,R); +else + [vbar,s,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R); +end +tEigs = cputime - tEigs; + + +%afficheTexte(sprintf('\n\nTemps ecoule pendant eigs : %g',tEigs),dataGraphCut.verbose,2); +%afficheTexte(sprintf('\nTemps ecoule pendant chol(Ubar''*Ubar) : %g',t_chol_Ubar),dataGraphCut.verbose); +if convergence~=0 + afficheTexte(sprintf(' (Non-convergence)'),dataGraphCut.verbose); +end + + +%disp(sprintf('nnz(P) : %f\n',nnz(P))); +%disp(sprintf('nnz(Ubar) : %f\n',nnz(Ubar))); +%disp(sprintf('nnz(R) : %f\n',nnz(R))); +%disp(sprintf('nnz(global) : %f\n',nnz(P) + 4 * nnz(Ubar) + 4*nnz(R))); + + + +s = real(diag(s)); +[x,y] = sort(-s); +s = -x; +vbar = vbar(:,y); + + +v = spdiags(Dinvsqrt,0,n,n) * vbar; + +for i=1:size(v,2) + %v(:,i) = v(:,i) / max(abs(v(:,i))); + v(:,i) = (v(:,i) / norm(v(:,i)) )*norm(ones(n,1)); + if v(1,i)~=0 + v(:,i) = - v(:,i) * sign(v(1,i)); + end +end + +% % nouvelle idee : factorisation de Cholesky +% t_chol_Ubar = cputime; +% R = chol(Ubar' * Ubar); +% t_chol_Ubar = cputime - t_chol_Ubar; +% R = sparsifyc(R,dataGraphCut.valeurMin); +% tEigs = cputime; +% if options.issym +% [vbar,s,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,R); +% else +% [vbar,s,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R); +% end +% tEigs = cputime - tEigs; + + +% % compute H = (Ubar'*Ubar)^(-1) +% t_inv_H = cputime; +% H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); %changer +% t_inv_H = cputime - t_inv_H; +% H = sparsifyc(H,dataGraphCut.valeurMin); +% tEigs = cputime; +% if options.issym & max(max(H-H')) < 1e-10 +% [vbar,s,convergence] = tim_eigs(@mex_projection_inv_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,triu(H)); +% else +% [vbar,s,convergence] = tim_eigs(@mex_projection_inv,n,nbEigenValues,'lm',options,P,Ubar,H); +% end +% tEigs = cputime - tEigs; + + + +% % idee de mon rapport... semble pas marcher car R = qr(Ubar,0) est plus +% % lent que H = inv(sparsifyc(Ubar' * Ubar,dataGraphCut.valeurMin)); +% t_qr_Ubar = cputime; +% R = qr(Ubar,0); +% t_qr_Ubar = cputime - t_qr_Ubar; +% R = sparsifyc(R,dataGraphCut.valeurMin); +% tEigs2 = cputime; +% if options.issym +% [vbar2,s2,convergence] = tim_eigs(@mex_projection_QR_symmetric,n,nbEigenValues,'lm',options,triu(P),Ubar,R); +% else +% [vbar2,s2,convergence] = tim_eigs(@mex_projection_QR,n,nbEigenValues,'lm',options,P,Ubar,R); +% end +% tEigs2 = cputime - tEigs2; + + + +% idee de Jianbo... semble pas marcher car on a besoin de prendre k maximal +% dans [A,S,B] = svds(Ubar,k); +% +% [A,S,B] = svds(Ubar,300); +% A = sparsifyc(A,dataGraphCut.valeurMin); +% tEigs = cputime; +% [vbar,s,convergence] = tim_eigs(@mex_projection_svd,n,nbEigenValues,'lm',options,P,A); + +% afficheTexte(sprintf('\ninv(H) : %g',t_inv_H),dataGraphCut.verbose); +% afficheTexte(sprintf('\n\nTemps ecoule pendant eigs : %g',tEigs2),dataGraphCut.verbose,2); +% afficheTexte(sprintf('\nqr(Ubar) : %g',t_qr_Ubar),dataGraphCut.verbose); +% disp(sprintf('nnz(H) : %f\n',nnz(H))); +% disp(sprintf('nnz(global) : %f\n',nnz(P) + 4 * nnz(Ubar) + 2*nnz(H))); diff --git a/SD-VBS/common/toolbox/MultiNcut/read_data.m b/SD-VBS/common/toolbox/MultiNcut/read_data.m new file mode 100755 index 0000000..c0929c0 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/read_data.m @@ -0,0 +1,13 @@ + +%fnames = dir('/home/jshi/Results_DLIB/SegLabl*.mat'); + +fnames = dir('/data/jshi/DLIB/Results/Results_DLIB/SegLabl*.mat'); + +for j=1:length(fnames), + cm = sprintf('load /data/jshi/DLIB/Results/Results_DLIB/%s',fnames(j).name); + disp(cm);eval(cm); +figure(1);imagesc(I); colormap(gray); axis image; +figure(2); imagesc(SegLabel); axis image; + +pause; +end diff --git a/SD-VBS/common/toolbox/MultiNcut/readimage.m b/SD-VBS/common/toolbox/MultiNcut/readimage.m new file mode 100755 index 0000000..e45fa19 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/readimage.m @@ -0,0 +1,15 @@ +function I = readimage(fn,maxSize); + +Io = imread(fn); +[nr,nc,nb] = size(Io); + +if nb>1, + I = rgb2gray(Io); +else + I= Io; +end + +%maxSize = 400; +if max(nr,nc) > maxSize, + I = imresize(I,maxSize/max(nr,nc)); +end diff --git a/SD-VBS/common/toolbox/MultiNcut/run_script.m b/SD-VBS/common/toolbox/MultiNcut/run_script.m new file mode 100755 index 0000000..3d2e8c6 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/run_script.m @@ -0,0 +1,60 @@ +%% set path for the MNcut code + +if 1, +% MNcutDir = '/home/jshi/Matlab/Toolbox/MultiNcut'; + MNcutDir = 'C:\qihuizhu\Checkout\Human\Source\MultiNcut_new\MultiNcut'; + path(path,MNcutDir); +end + +%% set the image input and output dir. +% imagedir = '/data/jshi/DLIB/image.cd'; +% imagedir = 'C:\qihuizhu\Checkout\Human\Source\Data\test'; +imagedir = 'C:\qihuizhu\Checkout\Human\Data\Current\baby_case5'; +% imageformat = 'ppm'; +imageformat = 'tif'; + +% OutputDir = '/home/jshi/Results_DLIB'; +% OutputDir = 'C:\qihuizhu\Checkout\Human\Source\Data\test'; +OutputDir = 'C:\qihuizhu\Checkout\Human\Result\Segmentation\MultiNcut_new_03.07'; + +a = dir(OutputDir); +if (length(a) == 0), + cm = sprintf('mkdir %s',OutputDir); + disp(cm); eval(cm); +end + +files = dir(sprintf('%s/*.%s',imagedir,imageformat)); + +%% image size definition +imageSize = 400; + +% for id =11:200, +for id = 1:length(files) + %for id = 19:19, + I=readimage(sprintf('%s/%s',imagedir,files(id).name),imageSize); + + num_segs = [10, 20]; + + tic + [SegLabel,eigenVectors,eigenValues]= MNcut(I,num_segs); + toc + + for j=1:size(SegLabel,3), + [gx,gy] = gradient(SegLabel(:,:,j)); + bw = (abs(gx)>0.1) + (abs(gy) > 0.1); + + figure(1);clf; J1=showmask(double(I),bw); imagesc(J1);axis image; + cm = sprintf('print -djpeg %s/file%.4d-%.2d.jpg',OutputDir,id,num_segs(j)); disp(cm);eval(cm); + + + figure(10);imagesc(SegLabel(:,:,j));axis image; + cm = sprintf('print -djpeg %s/Seg%.4d-%.2d.jpg',OutputDir,id,num_segs(j));disp(cm);eval(cm); + +% pause; + end + + fname = files(id).name; + %cm = sprintf('save %s/SegLabl%.4d.mat I SegLabel fname',OutputDir,id); disp(cm); eval(cm); + %cm = sprintf('save %s/SegEig%.4d.mat eigenVectors eigenValues',OutputDir,id);disp(cm); eval(cm); + +end diff --git a/SD-VBS/common/toolbox/MultiNcut/showmask.m b/SD-VBS/common/toolbox/MultiNcut/showmask.m new file mode 100755 index 0000000..6fd1142 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/showmask.m @@ -0,0 +1,65 @@ +% function RGB=showmask(V,M,hue); +% Input: +% V = image +% M = nonnegative mask +% hue = a number in [0,1], red,yellow,green,cyan,...,red +% a char, 'r','g','b','y','c','m' +% or a matrix of the same size of image +% eg. hue = mask1 * 0.7 + mask2 * 1; +% +% Output: +% RGB = an RGB image with V as shades and M as saturated hues +% If no output is required, this image is displayed. + +% Stella X. YU, 2000. Based on Jianbo Shi's version. + +function RGB=showmask(V,M,hue); + +if nargin<3 | isempty(hue), + hue = 0; +end +if ischar(hue), + switch hue, + case 'r', hue = 1.0; + case 'g', hue = 0.3; + case 'b', hue = 0.7; + case 'y', hue = 0.15; + case 'c', hue = 0.55; + case 'm', hue = 0.85; + end +end + + +V=V-min(V(:)); +V=V/max(V(:)); +V=.25+0.75*V; %brighten things up a bit + +M = double(M); +M = M-min(M(:)); +M = M/max(M(:)); + +H = hue+zeros(size(V)); +S = M; +RGB = hsv2rgb(H,S,V); + +if nargout>0, + return; +end + +hold off; image(RGB); axis('image'); +s = cell(1,2); +if isempty(inputname(1)), + s{1} = 'image'; +else + s{1} = inputname(1); +end +if isempty(inputname(2)), + s{2} = 'mask'; +else + s{2} = inputname(2); +end +title(sprintf('%s and colored %s',s{1},s{2})); + +if nargout==0, + clear RGB; +end diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c new file mode 100755 index 0000000..82fca98 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c @@ -0,0 +1,232 @@ +/*================================================================= +* syntax: SPMX = SPARSIFY(MX, THRES) +* +* SPARSIFY - sparsify the input matrix, i.e. ignore the values +* of the matrix which are below a threshold +* +* Input: - MX: m-by-n matrix (sparse or full) +* - THRES: threshold value (double) +* +* Output: - SPMX: m-by-n sparse matrix only with values +* whose absolut value is above the given threshold +* +* Written by Mirko Visontai (10/24/2003) +*=================================================================*/ + + +#include +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + /* Declare variable */ + int i,m,n,nzmax,newnnz,col,processed,passed; + int starting_row_index, current_row_index, stopping_row_index; + double *in_pr,*in_pi,*out_pr,*out_pi; + int *in_ir,*in_jc,*out_ir,*out_jc; + double thres; + + /* Check for proper number of input and output arguments */ + if ((nlhs != 1) || (nrhs != 2)){ + mexErrMsgTxt("usage: SPMX = SPARSIFY(MX, THRES)."); + } + /* if matrix is complex threshold the norm of the numbers */ + if (mxIsComplex(prhs[0])){ + /* Check data type of input argument */ + if (mxIsSparse(prhs[0])){ + + /* read input */ + in_pr = mxGetPr(prhs[0]); + in_pi = mxGetPi(prhs[0]); + in_ir = mxGetIr(prhs[0]); + in_jc = mxGetJc(prhs[0]); + nzmax = mxGetNzmax(prhs[0]); + m = mxGetM(prhs[0]); + n = mxGetN(prhs[0]); + thres = mxGetScalar(prhs[1]); + + /* Count new nonzeros */ + newnnz=0; + for(i=0; ithres) {newnnz++;} + } + + if (newnnz>0){ + /* create output */ + plhs[0] = mxCreateSparse(m,n,newnnz,mxCOMPLEX); + if (plhs[0]==NULL) + mexErrMsgTxt("Could not allocate enough memory!\n"); + out_pr = mxGetPr(plhs[0]); + out_pi = mxGetPr(plhs[0]); + out_ir = mxGetIr(plhs[0]); + out_jc = mxGetJc(plhs[0]); + passed = 0; + out_jc[0] = 0; + for (col=0; col thres){ + + out_pr[passed]=in_pr[current_row_index]; + out_pi[passed]=in_pi[current_row_index]; + out_ir[passed]=in_ir[current_row_index]; + out_jc[col+1] = out_jc[col+1]+1; + passed++; + } + } + } + } + } + else{ + plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX); + } + } + else{ /* for full matrices */ + /* read input */ + in_pr = mxGetPr(prhs[0]); + in_pi = mxGetPr(prhs[0]); + m = mxGetM(prhs[0]); + n = mxGetN(prhs[0]); + thres = mxGetScalar(prhs[1]); + + /* Count new nonzeros */ + newnnz=0; + for(i=0; ithres) {newnnz++;} + } + + if (newnnz>0){ + /* create output */ + plhs[0] = mxCreateSparse(m,n,newnnz,mxCOMPLEX); + if (plhs[0]==NULL) + mexErrMsgTxt("Could not allocate enough memory!\n"); + out_pr = mxGetPr(plhs[0]); + out_pi = mxGetPi(plhs[0]); + out_ir = mxGetIr(plhs[0]); + out_jc = mxGetJc(plhs[0]); + passed = 0; + out_jc[0] = 0; + + for (col=0; col thres){ + + out_pr[passed]=in_pr[current_row_index+m*col]; + out_ir[passed]=current_row_index; + out_jc[col+1] = out_jc[col+1]+1; + passed++; + } + } + } + } + else{ + plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX); + } + } + } + else { + /* Check data type of input argument */ + if (mxIsSparse(prhs[0])){ + + /* read input */ + in_pr = mxGetPr(prhs[0]); + in_ir = mxGetIr(prhs[0]); + in_jc = mxGetJc(prhs[0]); + nzmax = mxGetNzmax(prhs[0]); + n = mxGetN(prhs[0]); + m = mxGetM(prhs[0]); + thres = mxGetScalar(prhs[1]); + + /* Count new nonzeros */ + newnnz=0; + for(i=0; ithres) {newnnz++;} + } + + if (newnnz>0){ + /* create output */ + plhs[0] = mxCreateSparse(m,n,newnnz,mxREAL); + if (plhs[0]==NULL) + mexErrMsgTxt("Could not allocate enough memory!\n"); + out_pr = mxGetPr(plhs[0]); + out_ir = mxGetIr(plhs[0]); + out_jc = mxGetJc(plhs[0]); + passed = 0; + out_jc[0] = 0; + for (col=0; colthres){ + out_pr[passed]=in_pr[current_row_index]; + out_ir[passed]=in_ir[current_row_index]; + out_jc[col+1] = out_jc[col+1]+1; + passed++; + } + } + } + } + } + else{ + plhs[0] = mxCreateSparse(m,n,0,mxREAL); + } + } + else{ /* for full matrices */ + /* read input */ + in_pr = mxGetPr(prhs[0]); + n = mxGetN(prhs[0]); + m = mxGetM(prhs[0]); + thres = mxGetScalar(prhs[1]); + + /* Count new nonzeros */ + newnnz=0; + for(i=0; ithres) {newnnz++;} + } + + if (newnnz>0){ + /* create output */ + plhs[0] = mxCreateSparse(m,n,newnnz,mxREAL); + if (plhs[0]==NULL) + mexErrMsgTxt("Could not allocate enough memory!\n"); + out_pr = mxGetPr(plhs[0]); + out_ir = mxGetIr(plhs[0]); + out_jc = mxGetJc(plhs[0]); + passed = 0; + out_jc[0] = 0; + + for (col=0; colthres){ + out_pr[passed]=in_pr[current_row_index+m*col]; + out_ir[passed]=current_row_index; + out_jc[col+1] = out_jc[col+1]+1; + passed++; + } + } + } + } + else{ + plhs[0] = mxCreateSparse(m,n,0,mxREAL); + } + } + } +} + diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll new file mode 100755 index 0000000..cf832a6 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.dll differ diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64 b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64 new file mode 100755 index 0000000..2f5ed26 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexa64 differ diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx new file mode 100755 index 0000000..0ba41b6 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexglx differ diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac new file mode 100755 index 0000000..caa2306 Binary files /dev/null and b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.mexmac differ diff --git a/SD-VBS/common/toolbox/MultiNcut/spmtimesd.c b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.c new file mode 100755 index 0000000..a98dc0a --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/spmtimesd.c @@ -0,0 +1,141 @@ +/*================================================================ +* spmtimesd.c +* This routine computes a sparse matrix times a diagonal matrix +* whose diagonal entries are stored in a full vector. +* +* Examples: +* spmtimesd(m,d,[]) = diag(d) * m, +* spmtimesd(m,[],d) = m * diag(d) +* spmtimesd(m,d1,d2) = diag(d1) * m * diag(d2) +* m could be complex, but d is assumed real +* +* Stella X. Yu's first MEX function, Nov 9, 2001. + +% test sequence: + +m = 1000; +n = 2000; +a=sparse(rand(m,n)); +d1 = rand(m,1); +d2 = rand(n,1); +tic; b=spmtimesd(a,d1,d2); toc +tic; bb = spdiags(d1,0,m,m) * a * spdiags(d2,0,n,n); toc +e = (bb-b); +max(abs(e(:))) + +*=================================================================*/ + +# include "mex.h" + +void mexFunction( + int nargout, + mxArray *out[], + int nargin, + const mxArray *in[] +) +{ + /* declare variables */ + int i, j, k, m, n, nzmax, cmplx, xm, yn; + int *pir, *pjc, *qir, *qjc; + double *x, *y, *pr, *pi, *qr, *qi; + + /* check argument */ + if (nargin != 3) { + mexErrMsgTxt("Three input arguments required"); + } + if (nargout>1) { + mexErrMsgTxt("Too many output arguments."); + } + if (!(mxIsSparse(in[0]))) { + mexErrMsgTxt("Input argument #1 must be of type sparse"); + } + if ( mxIsSparse(in[1]) || mxIsSparse(in[2]) ) { + mexErrMsgTxt("Input argument #2 & #3 must be of type full"); + } + + /* computation starts */ + m = mxGetM(in[0]); + n = mxGetN(in[0]); + pr = mxGetPr(in[0]); + pi = mxGetPi(in[0]); + pir = mxGetIr(in[0]); + pjc = mxGetJc(in[0]); + + i = mxGetM(in[1]); + j = mxGetN(in[1]); + xm = ((i>j)? i: j); + + i = mxGetM(in[2]); + j = mxGetN(in[2]); + yn = ((i>j)? i: j); + + if ( xm>0 && xm != m) { + mexErrMsgTxt("Row multiplication dimension mismatch."); + } + if ( yn>0 && yn != n) { + mexErrMsgTxt("Column multiplication dimension mismatch."); + } + + + nzmax = mxGetNzmax(in[0]); + cmplx = (pi==NULL ? 0 : 1); + out[0] = mxCreateSparse(m,n,nzmax,cmplx); + if (out[0]==NULL) { + mexErrMsgTxt("Not enough space for the output matrix."); + } + + qr = mxGetPr(out[0]); + qi = mxGetPi(out[0]); + qir = mxGetIr(out[0]); + qjc = mxGetJc(out[0]); + + /* left multiplication */ + x = mxGetPr(in[1]); + if (yn==0) { + for (j=0; j 3) + error('Too many output arguments.') +end + +% Process inputs and do error-checking +if isa(varargin{1},'double') + A = varargin{1}; + Amatrix = 1; +else + A = fcnchk(varargin{1}); + Amatrix = 0; +end + +isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma) +if Amatrix + isrealprob = isreal(A); +end + +issymA = 0; +if Amatrix + issymA = isequal(A,A'); +end + +if Amatrix + [m,n] = size(A); + if (m ~= n) + error('A must be a square matrix or a function which computes A*x.') + end +else + n = varargin{2}; + nstr = 'Size of problem, ''n'', must be a positive integer.'; + if ~isequal(size(n),[1,1]) | ~isreal(n) + error(nstr) + end + if (round(n) ~= n) + warning('MATLAB:eigs:NonIntegerSize',['%s\n ' ... + 'Rounding input size.'],nstr) + n = round(n); + end + if issparse(n) + n = full(n); + end +end + +Bnotthere = 0; +Bstr = sprintf(['Generalized matrix B must be the same size as A and' ... + ' either a symmetric positive (semi-)definite matrix or' ... + ' its Cholesky factor.']); +if (nargin < (3-Amatrix-Bnotthere)) + B = []; + Bnotthere = 1; +else + Bk = varargin{3-Amatrix-Bnotthere}; + if isempty(Bk) % allow eigs(A,[],k,sigma,opts); + B = Bk; + else + if isequal(size(Bk),[1,1]) & (n ~= 1) + B = []; + k = Bk; + Bnotthere = 1; + else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ... + B = Bk; + if ~isa(B,'double') | ~isequal(size(B),[n,n]) + error(Bstr) + end + isrealprob = isrealprob & isreal(B); + end + end +end + +if Amatrix & ((nargin - ~Bnotthere)>4) + error('Too many inputs.') +end + +if (nargin < (4-Amatrix-Bnotthere)) + k = min(n,6); +else + k = varargin{4-Amatrix-Bnotthere}; +end + +kstr = ['Number of eigenvalues requested, k, must be a' ... + ' positive integer <= n.']; +if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n) + error(kstr) +end +if issparse(k) + k = full(k); +end +if (round(k) ~= k) + warning('MATLAB:eigs:NonIntegerEigQty',['%s\n ' ... + 'Rounding number of eigenvalues.'],kstr) + k = round(k); +end + +whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string.']); +if (nargin < (5-Amatrix-Bnotthere)) + % default: eigs('LM') => ARPACK which='LM', sigma=0 + eigs_sigma = 'LM'; + whch = 'LM'; + sigma = 0; +else + eigs_sigma = varargin{5-Amatrix-Bnotthere}; + if isstr(eigs_sigma) + % eigs(string) => ARPACK which=string, sigma=0 + if ~isequal(size(eigs_sigma),[1,2]) + whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ... + 'LM','SM','LA','SA','BE')]; + whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ... + 'LM','SM','LR','SR','LI','SI')]; + error(whchstr) + end + eigs_sigma = upper(eigs_sigma); + if isequal(eigs_sigma,'SM') + % eigs('SM') => ARPACK which='LM', sigma=0 + whch = 'LM'; + else + % eigs(string), where string~='SM' => ARPACK which=string, sigma=0 + whch = eigs_sigma; + end + sigma = 0; + else + % eigs(scalar) => ARPACK which='LM', sigma=scalar + if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1]) + error('Eigenvalue shift sigma must be a scalar.') + end + sigma = eigs_sigma; + if issparse(sigma) + sigma = full(sigma); + end + isrealprob = isrealprob & isreal(sigma); + whch = 'LM'; + end +end + +tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS) +maxit = []; +p = []; +info = int32(0); % use a random starting vector +display = 1; +cholB = 0; + +if (nargin >= (6-Amatrix-Bnotthere)) + opts = varargin{6-Amatrix-Bnotthere}; + if ~isa(opts,'struct') + error('Options argument must be a structure.') + end + + if isfield(opts,'issym') & ~Amatrix + issymA = opts.issym; + if (issymA ~= 0) & (issymA ~= 1) + error('opts.issym must be 0 or 1.') + end + end + + if isfield(opts,'isreal') & ~Amatrix + if (opts.isreal ~= 0) & (opts.isreal ~= 1) + error('opts.isreal must be 0 or 1.') + end + isrealprob = isrealprob & opts.isreal; + end + + if ~isempty(B) & (isfield(opts,'cholB') | isfield(opts,'permB')) + if isfield(opts,'cholB') + cholB = opts.cholB; + if (cholB ~= 0) & (cholB ~= 1) + error('opts.cholB must be 0 or 1.') + end + if isfield(opts,'permB') + if issparse(B) & cholB + permB = opts.permB; + if ~isequal(sort(permB),(1:n)) & ... + ~isequal(sort(permB),(1:n)') + error('opts.permB must be a permutation of 1:n.') + end + else + warning('MATLAB:eigs:IgnoredOptionPermB', ... + ['Ignoring opts.permB since B is not its sparse' ... + ' Cholesky factor.']) + end + else + permB = 1:n; + end + end + end + + if isfield(opts,'tol') + if ~isequal(size(opts.tol),[1,1]) | ~isreal(opts.tol) | (opts.tol<=0) + error(['Convergence tolerance opts.tol must be a strictly' ... + ' positive real scalar.']) + else + tol = full(opts.tol); + end + end + + if isfield(opts,'p') + p = opts.p; + pstr = ['Number of basis vectors opts.p must be a positive' ... + ' integer <= n.']; + if ~isequal(size(p),[1,1]) | ~isreal(p) | (p<=0) | (p>n) + error(pstr) + end + if issparse(p) + p = full(p); + end + if (round(p) ~= p) + warning('MATLAB:eigs:NonIntegerVecQty',['%s\n ' ... + 'Rounding number of basis vectors.'],pstr) + p = round(p); + end + end + + if isfield(opts,'maxit') + maxit = opts.maxit; + str = ['Maximum number of iterations opts.maxit must be' ... + ' a positive integer.']; + if ~isequal(size(maxit),[1,1]) | ~isreal(maxit) | (maxit<=0) + error(str) + end + if issparse(maxit) + maxit = full(maxit); + end + if (round(maxit) ~= maxit) + warning('MATLAB:eigs:NonIntegerIterationQty',['%s\n ' ... + 'Rounding number of iterations.'],str) + maxit = round(maxit); + end + end + + if isfield(opts,'v0') + if ~isequal(size(opts.v0),[n,1]) + error('Start vector opts.v0 must be n-by-1.') + end + if isrealprob + if ~isreal(opts.v0) + error(['Start vector opts.v0 must be real for real problems.']) + end + resid = full(opts.v0); + else + resid(1:2:(2*n-1),1) = full(real(opts.v0)); + resid(2:2:2*n,1) = full(imag(opts.v0)); + end + info = int32(1); % use resid as starting vector + end + + if isfield(opts,'disp') + display = opts.disp; + dispstr = 'Diagnostic level opts.disp must be an integer.'; + if (~isequal(size(display),[1,1])) | (~isreal(display)) | (display<0) + error(dispstr) + end + if (round(display) ~= display) + warning('MATLAB:eigs:NonIntegerDiagnosticLevel', ... + '%s\n Rounding diagnostic level.',dispstr) + display = round(display); + end + end + + if isfield(opts,'cheb') + warning('MATLAB:eigs:ObsoleteOptionCheb', ... + ['Ignoring polynomial acceleration opts.cheb' ... + ' (no longer an option).']); + end + if isfield(opts,'stagtol') + warning('MATLAB:eigs:ObsoleteOptionStagtol', ... + ['Ignoring stagnation tolerance opts.stagtol' ... + ' (no longer an option).']); + end + +end + +% Now we know issymA, isrealprob, cholB, and permB + +if isempty(p) + if isrealprob & ~issymA + p = min(max(2*k+1,20),n); + else + p = min(max(2*k,20),n); + end +end +if isempty(maxit) + maxit = max(300,ceil(2*n/max(p,1))); +end +if (info == int32(0)) + if isrealprob + resid = zeros(n,1); + else + resid = zeros(2*n,1); + end +end + +if ~isempty(B) % B must be symmetric (Hermitian) positive (semi-)definite + if cholB + if ~isequal(triu(B),B) + error(Bstr) + end + else + if ~isequal(B,B') + error(Bstr) + end + end +end + +useeig = 0; +if isrealprob & issymA + knstr = sprintf(['For real symmetric problems, must have number' ... + ' of eigenvalues k < n.\n']); +else + knstr = sprintf(['For nonsymmetric and complex problems, must have' ... + ' number of eigenvalues k < n-1.\n']); +end +if isempty(B) + knstr = [knstr sprintf(['Using eig(full(A)) instead.'])]; +else + knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])]; +end +if (k == 0) + useeig = 1; +end +if isrealprob & issymA + if (k > n-1) + if (n >= 6) + warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ... + '%s',knstr) + end + useeig = 1; + end +else + if (k > n-2) + if (n >= 7) + warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ... + '%s',knstr) + end + useeig = 1; + end +end + +if isrealprob & issymA + if ~isreal(sigma) + error(['For real symmetric problems, eigenvalue shift sigma must' ... + ' be real.']) + end +else + if ~isrealprob & issymA & ~isreal(sigma) + warning('MATLAB:eigs:ComplexShiftForRealProblem', ... + ['Complex eigenvalue shift sigma on a Hermitian problem' ... + ' (all real eigenvalues).']) + end +end + +if isrealprob & issymA + if strcmp(whch,'LR') + whch = 'LA'; + warning('MATLAB:eigs:SigmaChangedToLA', ... + ['For real symmetric problems, sigma value ''LR''' ... + ' (Largest Real) is now ''LA'' (Largest Algebraic).']) + end + if strcmp(whch,'SR') + whch = 'SA'; + warning('MATLAB:eigs:SigmaChangedToSA', ... + ['For real symmetric problems, sigma value ''SR''' ... + ' (Smallest Real) is now ''SA'' (Smallest Algebraic).']) + end + if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'}) + whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ... + 'LM','SM','LA','SA','BE')]; + error(whchstr) + end +else + if strcmp(whch,'BE') + warning('MATLAB:eigs:SigmaChangedToLM', ... + ['Sigma value ''BE'' is now only available for real' ... + ' symmetric problems. Computing ''LM'' eigenvalues instead.']) + whch = 'LM'; + end + if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'}) + whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ... + 'LM','SM','LR','SR','LI','SI')]; + error(whchstr) + end +end + +% Now have enough information to do early return on cases eigs does not handle +if useeig + if (nargout <= 1) + varargout{1} = eigs2(A,n,B,k,whch,sigma,cholB, ... + varargin{7-Amatrix-Bnotthere:end}); + else + [varargout{1},varargout{2}] = eigs2(A,n,B,k,whch,sigma,cholB, ... + varargin{7-Amatrix-Bnotthere:end}); + end + if (nargout == 3) + varargout{3} = 0; % flag indicates "convergence" + end + return +end + +if isrealprob & ~issymA + sigmar = real(sigma); + sigmai = imag(sigma); +end + +if isrealprob & issymA + if (p <= k) + error(['For real symmetric problems, must have number of' ... + ' basis vectors opts.p > k.']) + end +else + if (p <= k+1) + error(['For nonsymmetric and complex problems, must have number of' ... + ' basis vectors opts.p > k+1.']) + end +end + +if isequal(whch,'LM') & ~isequal(eigs_sigma,'LM') + % A*x = lambda*M*x, M symmetric (positive) semi-definite + % => OP = inv(A - sigma*M)*M and B = M + % => shift-and-invert mode + mode = 3; +elseif isempty(B) + % A*x = lambda*x + % => OP = A and B = I + mode = 1; +else % B is not empty + % Do not use mode=2. + % Use mode = 1 with OP = R'\(A*(R\x)) and B = I + % where R is B's upper triangular Cholesky factor: B = R'*R. + % Finally, V = R\V returns the actual generalized eigenvectors of A and B. + mode = 1; +end + +if cholB + pB = 0; + RB = B; + RBT = B'; +end + +if (mode == 3) & Amatrix % need lu(A-sigma*B) + if issparse(A) & (isempty(B) | issparse(B)) + if isempty(B) + AsB = A - sigma * speye(n); + else + if cholB + AsB = A - sigma * RBT * RB; + else + AsB = A - sigma * B; + end + end + [L,U,P,Q] = lu(AsB); + [perm,dummy] = find(Q); + else + if isempty(B) + AsB = A - sigma * eye(n); + else + if cholB + AsB = A - sigma * RBT * RB; + else + AsB = A - sigma * B; + end + end + [L,U,P] = lu(AsB); + end + dU = diag(U); + rcondestU = full(min(abs(dU)) / max(abs(dU))); + if (rcondestU < eps) + if isempty(B) + ds = sprintf(['(A-sigma*I) has small reciprocal condition' ... + ' estimate: %f\n'],rcondestU); + else + ds = sprintf(['(A-sigma*B) has small reciprocal condition' ... + ' estimate: %f\n'],rcondestU); + end + ds = [ds sprintf(['indicating that sigma is near an exact' ... + ' eigenvalue. The\nalgorithm may not converge unless' ... + ' you try a new value for sigma.\n'])]; + warning('MATLAB:eigs:SigmaNearExactEig',ds) + end +end + +if (mode == 1) & ~isempty(B) & ~cholB % need chol(B) + if issparse(B) + permB = symamd(B); + [RB,pB] = chol(B(permB,permB)); + else + [RB,pB] = chol(B); + end + if (pB == 0) + RBT = RB'; + else + error(Bstr) + end +end + +% Allocate outputs and ARPACK work variables +if isrealprob + if issymA % real and symmetric + prefix = 'ds'; + v = zeros(n,p); + ldv = int32(size(v,1)); + ipntr = int32(zeros(15,1)); + workd = zeros(n,3); + lworkl = p*(p+8); + workl = zeros(lworkl,1); + lworkl = int32(lworkl); + d = zeros(k,1); + else % real but not symmetric + prefix = 'dn'; + v = zeros(n,p); + ldv = int32(size(v,1)); + ipntr = int32(zeros(15,1)); + workd = zeros(n,3); + lworkl = 3*p*(p+2); + workl = zeros(lworkl,1); + lworkl = int32(lworkl); + workev = zeros(3*p,1); + d = zeros(k+1,1); + di = zeros(k+1,1); + end +else % complex + prefix = 'zn'; + zv = zeros(2*n*p,1); + ldv = int32(n); + ipntr = int32(zeros(15,1)); + workd = complex(zeros(n,3)); + zworkd = zeros(2*prod(size(workd)),1); + lworkl = 3*p^2+5*p; + workl = zeros(2*lworkl,1); + lworkl = int32(lworkl); + workev = zeros(2*2*p,1); + zd = zeros(2*(k+1),1); + rwork = zeros(p,1); +end + +ido = int32(0); % reverse communication parameter +if isempty(B) | (~isempty(B) & (mode == 1)) + bmat = 'I'; % standard eigenvalue problem +else + bmat = 'G'; % generalized eigenvalue problem +end +nev = int32(k); % number of eigenvalues requested +ncv = int32(p); % number of Lanczos vectors +iparam = int32(zeros(11,1)); +iparam([1 3 7]) = int32([1 maxit mode]); +select = int32(zeros(p,1)); + +cputms(1) = cputime - t0; % end timing pre-processing + +iter = 0; +ariter = 0; + + +%tim + + +indexArgumentsAfun = 7-Amatrix-Bnotthere:length(varargin); +nbArgumentsAfun = length(indexArgumentsAfun); +if nbArgumentsAfun >=1 + arguments_Afun = varargin{7-Amatrix-Bnotthere}; +end +if nbArgumentsAfun >=2 + arguments_Afun2 = varargin{7-Amatrix-Bnotthere+1}; +end +if nbArgumentsAfun >=3 + arguments_Afun3 = varargin{7-Amatrix-Bnotthere+2}; +end +%fin tim + + + +while (ido ~= 99) + + t0 = cputime; % start timing ARPACK calls **aupd + + if isrealprob + arpackc( [prefix 'aupd'], ido, ... + bmat, int32(n), whch, nev, tol, resid, ncv, ... + v, ldv, iparam, ipntr, workd, workl, lworkl, info); + else + zworkd(1:2:end-1) = real(workd); + zworkd(2:2:end) = imag(workd); + arpackc( 'znaupd', ido, ... + bmat, int32(n), whch, nev, tol, resid, ncv, zv, ... + ldv, iparam, ipntr, zworkd, workl, lworkl, ... + rwork, info ); + workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]); + end + + if (info < 0) + es = sprintf('Error with ARPACK routine %saupd: info = %d',... + prefix,full(info)); + error(es) + end + + cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd + t0 = cputime; % start timing MATLAB OP(X) + + % Compute which columns of workd ipntr references + + + + + + %[row,col1] = ind2sub([n,3],double(ipntr(1))); + %tim + row = mod(double(ipntr(1))-1,n) + 1 ; + col1 = floor((double(ipntr(1))-1)/n) + 1; + + + if (row ~= 1) + str = sprintf(['ipntr(1)=%d does not refer to the start of a' ... + ' column of the %d-by-3 array workd.'],full(ipntr(1)),n); + error(str) + end + + + + %[row,col2] = ind2sub([n,3],double(ipntr(2))); + %tim + row = mod(double(ipntr(2))-1,n) + 1 ; + col2 = floor((double(ipntr(2))-1)/n) + 1; + + + + if (row ~= 1) + str = sprintf(['ipntr(2)=%d does not refer to the start of a' ... + ' column of the %d-by-3 array workd.'],full(ipntr(2)),n); + error(str) + end + if ~isempty(B) & (mode == 3) & (ido == 1) + [row,col3] = ind2sub([n,3],double(ipntr(3))); + if (row ~= 1) + str = sprintf(['ipntr(3)=%d does not refer to the start of a' ... + ' column of the %d-by-3 array workd.'],full(ipntr(3)),n); + error(str) + end + end + + switch (ido) + case {-1,1} + if Amatrix + if (mode == 1) + if isempty(B) + % OP = A*x + workd(:,col2) = A * workd(:,col1); + else + % OP = R'\(A*(R\x)) + if issparse(B) + workd(permB,col2) = RB \ workd(:,col1); + workd(:,col2) = A * workd(:,col2); + workd(:,col2) = RBT \ workd(permB,col2); + else + workd(:,col2) = RBT \ (A * (RB \ workd(:,col1))); + end + end + elseif (mode == 3) + if isempty(B) + if issparse(A) + workd(perm,col2) = U \ (L \ (P * workd(:,col1))); + else + workd(:,col2) = U \ (L \ (P * workd(:,col1))); + end + else % B is not empty + if (ido == -1) + if cholB + workd(:,col2) = RBT * (RB * workd(:,col1)); + else + workd(:,col2) = B * workd(:,col1); + end + if issparse(A) & issparse(B) + workd(perm,col2) = U \ (L \ (P * workd(:,col1))); + else + workd(:,col2) = U \ (L \ (P * workd(:,col1))); + end + elseif (ido == 1) + if issparse(A) & issparse(B) + workd(perm,col2) = U \ (L \ (P * workd(:,col3))); + else + workd(:,col2) = U \ (L \ (P * workd(:,col3))); + end + end + end + else % mode is not 1 or 3 + error(['Unknown mode returned from ' prefix 'aupd.']) + end + else % A is not a matrix + if (mode == 1) + if isempty(B) + % OP = A*x + %workd(:,col2) = feval(A,workd(:,col1),varargin{7-Amatrix-Bnotthere:end}); + + + + + + nombre_A_times_X = nombre_A_times_X + 1; + + + + pause(0); %voir + + if nbArgumentsAfun == 1 + workd(:,col2) = feval(A,workd(:,col1),arguments_Afun); + %workd(:,col2) = max(workd(:,col2),0); + elseif nbArgumentsAfun == 2 + workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2); + elseif nbArgumentsAfun == 3 + workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2,arguments_Afun3); + else + workd(:,col2) = feval(A,workd(:,col1),varargin{indexArgumentsAfun}); + end + %workd(:,col2) = tim_w_times_x_c(workd(:,col1),arguments_Afun); %slower + + else + % OP = R'\(A*(R\x)) + if issparse(B) + workd(permB,col2) = RB \ workd(:,col1); + workd(:,col2) = feval(A,workd(:,col2),arguments_Afun); + workd(:,col2) = RBT \ workd(permB,col2); + + else + workd(:,col2) = RBT \ feval(A,(RB\workd(:,col1)),arguments_Afun); + end + end + elseif (mode == 3) + if isempty(B) + workd(:,col2) = feval(A,workd(:,col1), arguments_Afun); + else + if (ido == -1) + if cholB + workd(:,col2) = RBT * (RB * workd(:,col1)); + else + workd(:,col2) = B * workd(:,col1); + end + workd(:,col2) = feval(A,workd(:,col2), arguments_Afun); + elseif (ido == 1) + workd(:,col2) = feval(A,workd(:,col3), arguments_Afun); + end + end + else % mode is not 1 or 3 + error(['Unknown mode returned from ' prefix 'aupd.']) + end + end % if Amatrix + case 2 + if (mode == 3) + if cholB + workd(:,col2) = RBT * (RB * workd(:,col1)); + else + workd(:,col2) = B * workd(:,col1); + end + else + error(['Unknown mode returned from ' prefix 'aupd.']) + end + case 3 + % setting iparam(1) = ishift = 1 ensures this never happens + warning('MATLAB:eigs:WorklShiftsUnsupported', ... + ['EIGS does not yet support computing the shifts in workl.' ... + ' Setting reverse communication parameter to 99 and returning.']) + ido = int32(99); + case 99 + otherwise + error(['Unknown value of reverse communication parameter' ... + ' returned from ' prefix 'aupd.']) + end + + cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X) + + %tim + if nombreIterations ~= double(ipntr(15)) + nombreIterations = double(ipntr(15)); + end + + if display >= 1 && display <=2 + iter = double(ipntr(15)); + if (iter > ariter) & (ido ~= 99) + ariter = iter; + ds = sprintf(['Iteration %d: a few Ritz values of the' ... + ' %d-by-%d matrix:'],iter,p,p); + disp(ds) + if isrealprob + if issymA + dispvec = [workl(double(ipntr(6))+(0:p-1))]; + if isequal(whch,'BE') + % roughly k Large eigenvalues and k Small eigenvalues + disp(dispvec(max(end-2*k+1,1):end)) + else + % k eigenvalues + disp(dispvec(max(end-k+1,1):end)) + end + else + dispvec = [complex(workl(double(ipntr(6))+(0:p-1)), ... + workl(double(ipntr(7))+(0:p-1)))]; + % k+1 eigenvalues (keep complex conjugate pairs together) + disp(dispvec(max(end-k,1):end)) + end + else + dispvec = [complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ... + workl(2*double(ipntr(6))+(0:2:2*(p-1))))]; + disp(dispvec(max(end-k+1,1):end)) + end + end + end + +end % while (ido ~= 99) + +t0 = cputime; % start timing post-processing + +flag = 0; +if (info < 0) + es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,full(info)); + error(es) +else + if (nargout >= 2) + rvec = int32(1); % compute eigenvectors + else + rvec = int32(0); % do not compute eigenvectors + end + + if isrealprob + if issymA + arpackc( 'dseupd', rvec, 'A', select, ... + d, v, ldv, sigma, ... + bmat, int32(n), whch, nev, tol, resid, ncv, ... + v, ldv, iparam, ipntr, workd, workl, lworkl, info ); + if isequal(whch,'LM') | isequal(whch,'LA') + d = flipud(d); + if (rvec == 1) + v(:,1:k) = v(:,k:-1:1); + end + end + if ((isequal(whch,'SM') | isequal(whch,'SA')) & (rvec == 0)) + d = flipud(d); + end + else + arpackc( 'dneupd', rvec, 'A', select, ... + d, di, v, ldv, sigmar, sigmai, workev, ... + bmat, int32(n), whch, nev, tol, resid, ncv, ... + v, ldv, iparam, ipntr, workd, workl, lworkl, info ); + d = complex(d,di); + if rvec + d(k+1) = []; + else + zind = find(d == 0); + if isempty(zind) + d = d(k+1:-1:2); + else + d(max(zind)) = []; + d = flipud(d); + end + end + end + else + zsigma = [real(sigma); imag(sigma)]; + arpackc( 'zneupd', rvec, 'A', select, ... + zd, zv, ldv, zsigma, workev, ... + bmat, int32(n), whch, nev, tol, resid, ncv, zv, ... + ldv, iparam, ipntr, zworkd, workl, lworkl, ... + rwork, info ); + if issymA + d = zd(1:2:end-1); + else + d = complex(zd(1:2:end-1),zd(2:2:end)); + end + v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]); + end + + if (info ~= 0) + es = ['Error with ARPACK routine ' prefix 'eupd: ']; + switch double(info) + case 2 + ss = sum(select); + if (ss < k) + es = [es ... + ' The logical variable select was only set with ' int2str(ss) ... + ' 1''s instead of nconv=' int2str(double(iparam(5))) ... + ' (k=' int2str(k) ').' ... + ' Please contact the ARPACK authors at arpack@caam.rice.edu.']; + else + es = [es ... + 'The LAPACK reordering routine ' prefix(1) ... + 'trsen did not return all ' int2str(k) ' eigenvalues.'] + end + case 1 + es = [es ... + 'The Schur form could not be reordered by the LAPACK routine ' ... + prefix(1) 'trsen.' ... + ' Please contact the ARPACK authors at arpack@caam.rice.edu.']; + case -14 + es = [es prefix ... + 'aupd did not find any eigenvalues to sufficient accuracy.']; + otherwise + es = [es sprintf('info = %d',full(info))]; + end + error(es) + else + nconv = double(iparam(5)); + if (nconv == 0) + if (nargout < 3) + warning('MATLAB:eigs:NoEigsConverged', ... + 'None of the %d requested eigenvalues converged.',k) + else + flag = 1; + end + elseif (nconv < k) + if (nargout < 3) + warning('MATLAB:eigs:NotAllEigsConverged', ... + 'Only %d of the %d requested eigenvalues converged.',nconv,k) + else + flag = 1; + end + end + end % if (info ~= 0) +end % if (info < 0) + +if (issymA) | (~isrealprob) + if (nargout <= 1) + if isrealprob + varargout{1} = d; + else + varargout{1} = d(k:-1:1,1); + end + else + varargout{1} = v(:,1:k); + varargout{2} = diag(d(1:k,1)); + if (nargout >= 3) + varargout{3} = flag; + end + end +else + if (nargout <= 1) + varargout{1} = d; + else + cplxd = find(di ~= 0); + % complex conjugate pairs of eigenvalues occur together + cplxd = cplxd(1:2:end); + v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ... + complex(v(:,cplxd),-v(:,cplxd+1))]; + varargout{1} = v(:,1:k); + varargout{2} = diag(d); + if (nargout >= 3) + varargout{3} = flag; + end + end +end + +if (nargout >= 2) & (mode == 1) & ~isempty(B) + if issparse(B) + varargout{1}(permB,:) = RB \ varargout{1}; + else + varargout{1} = RB \ varargout{1}; + end +end + +cputms(4) = cputime-t0; % end timing post-processing + +cputms(5) = sum(cputms(1:4)); % total time + +if (display >= 2) %tim + if (mode == 1) + innerstr = sprintf(['Compute A*X:' ... + ' %f\n'],cputms(3)); + elseif (mode == 2) + innerstr = sprintf(['Compute A*X and solve B*X=Y for X:' ... + ' %f\n'],cputms(3)); + elseif (mode == 3) + if isempty(B) + innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ... + ' %f\n'],cputms(3)); + else + innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ... + ' %f\n'],cputms(3)); + end + end + if ((mode == 3) & (Amatrix)) + if isempty(B) + prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ... + ' %f\n'],cputms(1)); + else + prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ... + ' %f\n'],cputms(1)); + end + elseif ((mode == 2) & (~cholB)) + prepstr = sprintf(['Pre-processing, including chol(B):' ... + ' %f\n'],cputms(1)); + else + prepstr = sprintf(['Pre-processing:' ... + ' %f\n'],cputms(1)); + end + sstr = sprintf(['***********CPU Timing Results in seconds***********']); + ds = sprintf(['\n' sstr '\n' ... + prepstr ... + 'ARPACK''s %saupd: %f\n' ... + innerstr ... + 'Post-processing with ARPACK''s %seupd: %f\n' ... + '***************************************************\n' ... + 'Total: %f\n' ... + sstr '\n'], ... + prefix,cputms(2),prefix,cputms(4),cputms(5)); + disp(ds) + if nombre_A_times_X > 0 %tim + disp(sprintf('Number of iterations : %f\n',nombreIterations)); + disp(sprintf('Number of times A*X was computed : %f\n',nombre_A_times_X)); + disp(sprintf('Average time for A*X : %f\n',cputms(3)/nombre_A_times_X)); + end +end -- cgit v1.2.2