vb_back_threshold

 Calculate threshold for background selection from MRI image
  [threshold] = vb_back_threshold(B)
  [threshold] = vb_back_threshold(B,pmax)
  [threshold, step, Nhist, pdf] = vb_back_threshold(B,pmax,Nwid,Imax)
 --- Input
 B    : MRI image
 --- Optional input
 pmax : Probability that the voxcel is background (= 0.998)
        if (image intensity) < threshold 
 Nwid : Width of Gaussian filter to smooth histgram
 Imax : Max number of histgram bin (256)
 --- Output
 threshold: threshold for background selection
   Background histogram is modeled by Rayleigh distribution
 --- Optional Output
 Nhist : Histogram of MRI image intensity
 pdf   : Estimated background histogram
 step  : Intensity level step size 

 if max(B) <= Imax, step = 1
 else               step = round(max(B)/Imax);

 Made by M. Sato 2004-3-28
 M. Sato 2010-10-10
  Histogram is smoothed by Gaussian filter with width (Nwid=5)
    to find first local minimum point
  Three histgram distributions (Rayleigh, Gaussian, Exponential) are examined
  if (error of Rayleigh) < 0.1, Rayleigh is selected
  othewise, least error method is selected
 M. Sato 2010-10-18
  Histogram level setting is modified

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function     [th, step, y, pdf, err] = vb_back_threshold(B,pmax,Nwid,Imax)
0002 % Calculate threshold for background selection from MRI image
0003 %  [threshold] = vb_back_threshold(B)
0004 %  [threshold] = vb_back_threshold(B,pmax)
0005 %  [threshold, step, Nhist, pdf] = vb_back_threshold(B,pmax,Nwid,Imax)
0006 % --- Input
0007 % B    : MRI image
0008 % --- Optional input
0009 % pmax : Probability that the voxcel is background (= 0.998)
0010 %        if (image intensity) < threshold
0011 % Nwid : Width of Gaussian filter to smooth histgram
0012 % Imax : Max number of histgram bin (256)
0013 % --- Output
0014 % threshold: threshold for background selection
0015 %   Background histogram is modeled by Rayleigh distribution
0016 % --- Optional Output
0017 % Nhist : Histogram of MRI image intensity
0018 % pdf   : Estimated background histogram
0019 % step  : Intensity level step size
0020 %
0021 % if max(B) <= Imax, step = 1
0022 % else               step = round(max(B)/Imax);
0023 %
0024 % Made by M. Sato 2004-3-28
0025 % M. Sato 2010-10-10
0026 %  Histogram is smoothed by Gaussian filter with width (Nwid=5)
0027 %    to find first local minimum point
0028 %  Three histgram distributions (Rayleigh, Gaussian, Exponential) are examined
0029 %  if (error of Rayleigh) < 0.1, Rayleigh is selected
0030 %  othewise, least error method is selected
0031 % M. Sato 2010-10-18
0032 %  Histogram level setting is modified
0033 %
0034 % Copyright (C) 2011, ATR All Rights Reserved.
0035 % License : New BSD License(see VBMEG_LICENSE.txt)
0036 
0037 % --------------------------------------------------------
0038 % MRI構造画像 B から背景消去の為の閾値設定
0039 % pmax : 背景ノイズが従うRayleigh分布に含まれると判断する閾値確率
0040 %        この値が大きいほど閾値が大きくなる
0041 %
0042 % 背景ノイズが従うRayleigh分布の分散を推定し、
0043 % この分布に含まれる確率が pmax である閾値を求める
0044 
0045 if ~exist('pmax','var'), pmax = 0.998; end;
0046 if ~exist('Imax','var'), Imax = 256; end;
0047 if ~exist('Nwid','var'), Nwid = 5; end;
0048 
0049 % Rayleigh threshold = rmax * (variance)
0050 rmax = sqrt(2*abs(log(1-pmax)));
0051 % Gaussian threshold = xpeak + rgas * (variance)
0052 rgas = erfinv(pmax) * sqrt(2);
0053 % Exponential threshold = xpeak + abs(rexp/c);
0054 % pdf = a + b * exp( - c * x );
0055 rexp = log(1-pmax);
0056 
0057 % Max error to select Rayleigh distribution
0058 err_max = 0.1;
0059 
0060 % y : MRI image histogram
0061 [y,x,step] = vb_mri_histgram(B(:), Imax);
0062 Bmax = max(B(:));
0063 
0064 % Smoothing histgram and find first local minimum
0065 %   z : smoothing histgram
0066 [ix_min, ix_max, z] = vb_histmin_estimate(y,Nwid);
0067 x0 = x(ix_min);
0068 
0069 % Max point of histogram
0070 [ymax ,imax]= max(y);
0071 if imax == 1, imax = 3; end;
0072 
0073 % Background histogram estimate region
0074 kmax = max(ceil(imax*rmax), ix_min);
0075 ix = 1:kmax;
0076 
0077 % Variance of Rayleigh distribution
0078 [s1,pdf1,err1] = vb_rayleigh_estimate(y,x,kmax);
0079 % Rayleigh threshold,
0080 x1 = s1*rmax;
0081 
0082 % Center and variance of Gaussian
0083 [x2,s2,C,err2] = vb_gaussian_estimate(y(imax:kmax),x(imax:kmax));
0084 
0085 % pdf2  : Gaussian distribution
0086 pdf2 = C * exp(-0.5*((x(ix) - x2)/s2).^2);
0087 
0088 % threshold,
0089 x2  = x2 + s2*rgas;
0090 
0091 % parameter of Exponential distribution
0092 
0093 [a, b, c, err3] = vb_exponential_estimate(y(imax:kmax), x(imax:kmax));
0094 
0095 % pdf3 : Exponential distribution
0096 pdf3 = a + b * exp( c * x(ix) );
0097 % ∫ c * exp( c*x ) = exp( c*x )
0098 x3  = x(imax) + abs(rexp/c);
0099 
0100 if err1 < err_max
0101     th  = x1;
0102     pdf = pdf1;
0103     err = err1;
0104     str = 'Rayleigh';
0105 else
0106     
0107     [tmp,jmode] = min([err2, err3]);
0108     
0109     switch    jmode
0110     case    1
0111         th  = x2;
0112         pdf = pdf2;
0113         err = err2;
0114         str = 'Gaussian';
0115     case    2
0116         th  = x3;
0117         pdf = pdf3;
0118         err = err3;
0119         str = 'Exponential';
0120     end
0121 end
0122 
0123 fprintf('Max intensity = %g\n',Bmax)
0124 fprintf('Histgram step = %d\n',step)
0125 
0126 fprintf('Rayleigh    = %g (err=%g)\n',x1,err1)
0127 fprintf('Gaussian    = %g (err=%g)\n',x2,err2)
0128 fprintf('Exponential = %g (err=%g)\n\n',x3,err3)
0129 
0130 fprintf('%s distribution estimation\n',str)
0131 fprintf('Threshold   = %g (err=%g)\n',th,err)
0132 
0133 if nargout < 5, return; end;
0134 %%%%%%%%%%%%%%%%%%%%% Plot Histgram %%%%%%%%%%%%%%%%%%%%
0135 
0136 zmax = max(z);
0137 zmid = z(ix_max);
0138 
0139 %ix_th = find( x >= th);
0140 ix_th = round(th/step);
0141 
0142 % Plot histgram
0143 subplot(2,2,1)
0144 plot(y)
0145 hold on
0146 plot(z,'-r')
0147 
0148 plot(ix, pdf1,'-c')
0149 plot(ix, pdf2,'-m')
0150 plot(ix, pdf3,'-y')
0151 
0152 plot([ix_th ix_th], [0 0.5*ymax],'-r')
0153 
0154 ylim([0 ymax]);
0155 xlim([0 ix_min])
0156 title('Low intensity histogram')
0157 
0158 % Plot histgram
0159 subplot(2,2,2)
0160 plot(y)
0161 hold on
0162 plot(z,'-r')
0163 
0164 plot(ix, pdf1,'-c')
0165 
0166 plot([ix_th ix_th], [0 8*zmax],'-r')
0167 
0168 ylim([0 zmax]);
0169 xlim([0 ix_min])
0170 title('Low intensity histogram')
0171 
0172 % Plot histgram
0173 subplot(2,2,3)
0174 plot(y)
0175 hold on
0176 plot(z,'-r')
0177 
0178 plot(ix, pdf2,'-m')
0179 plot(ix, pdf3,'-y')
0180 
0181 plot([ix_th ix_th], [0 8*zmax],'-r')
0182 
0183 ylim([0 zmax]);
0184 xlim([0 ix_min])
0185 title('Low intensity histogram')
0186 
0187 % Plot histgram
0188 %figure
0189 subplot(2,2,4)
0190 plot(y)
0191 hold on
0192 plot(z,'-r')
0193 
0194 plot(ix, pdf1,'-c')
0195 plot(ix, pdf2,'-m')
0196 plot(ix, pdf3,'-y')
0197 
0198 plot([ix_th ix_th], [0 zmax],'-r')
0199 
0200 ylim([0 2*zmid]);
0201 xlim([0 1.5*ix_max])
0202 title('Intensity histogram')
0203 
0204 return
0205 %%%%%%%%%%%%%%%%%%%%%%%%  END %%%%%%%%%%%%%%%%%%%%
0206 xx  = x(1:kmax);
0207 
0208 % Plot histgram
0209 subplot(2,2,1)
0210 plot(x,y)
0211 hold on
0212 plot(x,z,'-r')
0213 
0214 plot(xx, pdf1,'-c')
0215 plot(xx, pdf2,'-y')
0216 plot(xx, pdf3,'-m')
0217 
0218 plot([th th], [0 0.5*ymax],'-r')
0219 
0220 ylim([0 max(y)]);
0221 xlim([0 x0])
0222 title('Low intensity histogram')
0223 
0224 % Plot histgram
0225 subplot(2,2,2)
0226 plot(x,y)
0227 hold on
0228 plot(x,z,'-r')
0229 
0230 plot(xx, pdf1,'-c')
0231 
0232 plot([th th], [0 8*zmax],'-r')
0233 
0234 ylim([0 max(z)]);
0235 xlim([0 x0])
0236 title('Low intensity histogram')
0237 
0238 % Plot histgram
0239 subplot(2,2,3)
0240 plot(x,y)
0241 hold on
0242 plot(x,z,'-r')
0243 
0244 plot(xx, pdf2,'-c')
0245 plot(xx, pdf3,'-y')
0246 
0247 plot([th th], [0 8*zmax],'-r')
0248 
0249 ylim([0 max(z)]);
0250 xlim([0 x0])
0251 title('Low intensity histogram')
0252 
0253 % Plot histgram
0254 %figure
0255 subplot(2,2,4)
0256 plot(x,y)
0257 hold on
0258 plot(x,z,'-r')
0259 
0260 plot(xx, pdf,'-y')
0261 plot([th th], [0 zmax],'-r')
0262 
0263 ylim([0 2*zmax]);
0264 xlim([0 1.5*x(ix_max)])
0265 title('Low intensity histogram')
0266 title('Intensity histogram')
0267 
0268 return
0269 %%%%%%%%%%%%%%%%%%%%%%  END  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0270 
0271 Bhigh = max(Nhist(j0+1:end));
0272 
0273 subplot(2,1,1)
0274 
0275 plot(level,Nhist);
0276 hold on
0277 %plot([level(imax) level(imax)], [0 Hmax*1.1], 'b-');
0278 plot([th th], [0 0.5*max(Nhist)], 'r-');
0279 plot(x0,rpdf,'-r')
0280 xlim([0 2*th]);
0281 
0282 title('Low intensity histogram of MRI (Rayleigh dist)')
0283 
0284 subplot(2,1,2)
0285 plot(level,Nhist);
0286 hold on
0287 plot([th th], [0 Nhist(j0)], 'r-');
0288 plot(x0,rpdf,'-r')
0289 xlim([th*0.8 Bmax]);
0290 ylim([0 Bhigh]);
0291 
0292 title('Intensity histogram of MRI (Rayleigh dist)')
0293 
0294 return
0295 
0296 % Variance of Rayleigh distribution (smoothed histgram)
0297 [s2,pdf2,err2] = vb_rayleigh_estimate(z,x,kmax);
0298 % Rayleigh threshold,
0299 x2 = s2*rmax;
0300 
0301 % Find histogram increasing point
0302 jx = find( x >= x1);
0303 j0 = jx(1);
0304 dh = y(j0+1:end) - y(j0:end-1);
0305 jj = find( dh > 0 );
0306 j0 = j0 + jj(1) - 1;
0307 xm = x(j0);