From d17b33131c14864bd1eae275f49a3f148e21cf29 Mon Sep 17 00:00:00 2001 From: Leo Chan Date: Thu, 22 Oct 2020 01:53:21 -0400 Subject: Squashed commit of the sb-vbs branch. Includes the SD-VBS benchmarks modified to: - Use libextra to loop as realtime jobs - Preallocate memory before starting their main computation - Accept input via stdin instead of via argc Does not include the SD-VBS matlab code. Fixes libextra execution in LITMUS^RT. --- 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=gaussia