diff --git a/functions/readDicomStack.m b/functions/readDicomStack.m
index 23c82f8..69b33a9 100755
--- a/functions/readDicomStack.m
+++ b/functions/readDicomStack.m
@@ -1,353 +1,348 @@
 function [img, imgInfo] = readDicomStack(dicomDir, opt)
 %
 dicomdict('set','dicom.dic');
 
 % kernelTypes = {'bone', 'lung', 'none'};
 if(~isempty(dicomDir) && dicomDir(end) ~= filesep)
     dicomDir = [dicomDir filesep];
 end
 
 
 if (nargin == 0)
     dicomDir = './';
     opt = 'max';
 else
     if (nargin == 1)
         opt = 'max';
     end
 end
 
 if class(dicomDir)=='char'
     filesList = dir([dicomDir '*.dcm']);
     if(isempty(filesList))
         filesList = dir(dicomDir);
         filesList = filesList(3:end);
     end
     tmp{size(filesList, 1)} = '';
     for ii = 1:size(filesList, 1)
         tmp{ii} = [dicomDir filesList(ii).name];
     end
     filesList = tmp;
 else
     filesList = dicomDir;
 end
 
 
 Headers{size(filesList, 2)} = [];
 sliceLoc = zeros(size(filesList, 1), 1);
 kernelType = zeros(size(sliceLoc));
 for ii = 1:size(filesList, 2)
     try
         Headers{ii}  = dicominfo([filesList{ii}]);
     catch
         sliceLoc(ii) = inf;
         continue;
     end
     try
         sliceLoc(ii) = Headers{ii}.ImagePositionPatient(3);
     catch
         sliceLoc(ii) = -inf;
     end
     if (isfield(Headers{ii}, 'ConvolutionKernel'))
         switch lower(Headers{ii}.ConvolutionKernel)
             case 'bone'
                 kernelType(ii) = 1;
             case 'lung'
                 kernelType(ii) = 2;
             case 'standard'
                 kernelType(ii) = 2;
             case 'soft'
                 kernelType(ii) = 2;
             otherwise
                 kernelType(ii) = 3;
         end
     else
         kernelType(ii) = 3;
     end
 end
 
 [kernelType, ind] = sort(kernelType);
 sliceLoc = sliceLoc(ind);
 Headers = Headers(:, ind);
 numStacks = 0;
 
 ii = 1;
 while (ii<size(Headers, 2))
     jj = ii;
     while(jj<=size(Headers, 2) && kernelType(jj) == kernelType(ii)), jj = jj+1; end
     tmp = sliceLoc(ii:jj-1);
     [sliceLoc(ii:jj-1), ind] = sort(tmp);
     tmp = Headers(:, ii:jj-1);
     Headers(:, ii:jj-1) = tmp(:, ind);
     ii = jj;
     numStacks = numStacks+1;
 end
 
 if (numStacks == 0)
     imgInfo = [];
     img     = [];
     return;
 end
 
 
 currentSlice = 1;
 img{numStacks} = [];
 for ii = 1:numStacks
     currentKernel = kernelType(currentSlice);
     while(sliceLoc(currentSlice)==-inf), currentSlice = currentSlice+1; end
     xMin = 0; xMax = 0; xSize = 0;
     yMin = 0; yMax = 0; ySize = 0;
     zMin = 0; zMax = 0; zSize = 0;
     PixelSpacing = [0, 0];
     SliceThickness = 0;
 
     if (isfield(Headers{currentSlice}, 'PixelSpacing'))
         PixelSpacing = Headers{currentSlice}.PixelSpacing;
     end
     %     if (isfield(Headers{currentSlice}, 'SliceThickness'))
     %         SliceThickness = Headers{currentSlice}.SliceThickness;
     %     end
     if (isfield(Headers{currentSlice}, 'ImagePositionPatient'))
         xMin = Headers{currentSlice}.ImagePositionPatient(1);
         yMin = Headers{currentSlice}.ImagePositionPatient(2);
         zMin = Headers{currentSlice}.ImagePositionPatient(3);
     end
     if (isfield(Headers{currentSlice}, 'Columns'))
         xSize = Headers{currentSlice}.Columns;
         xMax  = xMin + double(xSize - 1) * PixelSpacing(2);
     end
     if (isfield(Headers{currentSlice}, 'Rows'))
         ySize = Headers{currentSlice}.Rows;
         yMax  = yMin + double(ySize - 1) * PixelSpacing(1);
     end
 
     switch currentKernel
         case 1
             imgInfo(ii).ConvolutionKernel = 'Bone';
         case 2
             imgInfo(ii).ConvolutionKernel = 'Lung';
         otherwise
             imgInfo(ii).ConvolutionKernel = 'None';
     end
     while(currentSlice <= numel(kernelType) && currentKernel == kernelType(currentSlice))
         zSize = zSize+1;
         try
             tmp = dicomread(Headers{currentSlice});
             if (isfield(Headers{currentSlice}, 'RescaleSlope') && Headers{currentSlice}.RescaleSlope~=0)
                 Slope = Headers{currentSlice}.RescaleSlope;
             else
                 Slope = 1;
             end
             if (isfield(Headers{currentSlice}, 'RescaleIntercept'))
                 Intercept = Headers{currentSlice}.RescaleIntercept;
             else
                 Intercept = 0;
             end
         catch
             currentSlice = currentSlice + 1;
             continue;
         end
         img{ii}(:, :, zSize) = tmp*Slope + Intercept;
         if(SliceThickness == 0 && currentSlice<numel(kernelType) && kernelType(currentSlice+1) == currentKernel)
             SliceThickness = sliceLoc(currentSlice+1) - sliceLoc(currentSlice);
         else
             newLoc = zMin+(zSize-1)*SliceThickness;
             if(SliceThickness ~= 0 && newLoc ~= sliceLoc(currentSlice) && currentSlice<numel(kernelType) && kernelType(currentSlice+1) == currentKernel)
                 sliceLoc(currentSlice) = newLoc;
             end
         end
         zMax  = sliceLoc(currentSlice);
         currentSlice = currentSlice + 1;
     end
     imgInfo(ii).BoundingBox.xMin = xMin; imgInfo(ii).BoundingBox.xMax = xMax; imgInfo(ii).xSize = xSize;
     imgInfo(ii).BoundingBox.yMin = yMin; imgInfo(ii).BoundingBox.yMax = yMax; imgInfo(ii).ySize = ySize;
     imgInfo(ii).BoundingBox.zMin = zMin; imgInfo(ii).BoundingBox.zMax = zMax; imgInfo(ii).zSize = zSize;
     imgInfo(ii).SliceThickness = (zMax -zMin)/(zSize-1);
     imgInfo(ii).PixelSpacing = [(xMax - xMin)/double(xSize-1) (yMax - yMin)/double(ySize-1)];
 
 end
 
 ind = [];
-display(['TEST : length(img) = ' num2str(length(img))]);
-display(['TEST : numStacks = ' num2str(numStacks)]);
-display(['TEST : opt = ' num2str(opt)]);
 for ii = 1:numStacks
     if(isempty(img{ii}))
         ind = [ind, ii];
     end
 end
 
 
 img(:, ind) = [];
 imgInfo(:, ind) = [];
 
 switch opt
     case 'all'
     case 'first'
         img = img{1};
         imgInfo = imgInfo(1);
     case 'last'
         img = img{end};
         imgInfo = imgInfo(end);
     otherwise
         SZs = zeros(size(imgInfo, 2), 1);
         for tmp = 1:size(imgInfo, 2)
             SZs(tmp) = imgInfo(tmp).zSize;
         end
 
         tmp = max(SZs);
         ind = find(SZs == tmp);
         if(numel(ind)>1)
             for ii=1:numel(ind)
                 if(imgInfo(ind(ii)).ConvolutionKernel == 'Lung')
                     break;
                 end
             end
             if (ii>numel(ind))
                 ind = ind(1);
             else
                 ind = ind(ii);
             end
         end
-        display(['TEST : length(img) = ' num2str(length(img))]);
-        display(['TEST : ind = ' num2str(ind)]);
         imgInfo = imgInfo(ind);
 end
 
 
 % figure
 % for ii = 1:size(img{1}, 3)
 %     imshow(uint8(img{1}(:, :, ii)/4095*255));
 %     pause(0.1)
 % end
 
 
 
 % % % stackNo = 1;
 % % %
 % % %
 % % %
 % % % info = dicominfo([dicomDir filesList(1).name]);
 % % % tmp  = dicomread(info);
 % % % imgInfo(stackNo).BoundingBox.xMin = info.ImagePositionPatient(1);
 % % % imgInfo(stackNo).BoundingBox.yMin = info.ImagePositionPatient(2);
 % % % imgInfo(stackNo).BoundingBox.zMin = info.ImagePositionPatient(3);
 % % % imgInfo(stackNo).PixelSpacing     = info.PixelSpacing;
 % % % imgInfo(stackNo).SliceThickness   = 0;
 % % % imgInfo(stackNo).xSize = info.Columns;
 % % % imgInfo(stackNo).ySize = info.Rows;
 % % % imgInfo(stackNo).zSize = 1;
 % % % imgInfo(stackNo).BoundingBox.xMax = imgInfo(stackNo).BoundingBox.xMin + double(imgInfo(stackNo).xSize -1) * imgInfo(stackNo).PixelSpacing(1);
 % % % imgInfo(stackNo).BoundingBox.yMax = imgInfo(stackNo).BoundingBox.yMin + double(imgInfo(stackNo).ySize -1) * imgInfo(stackNo).PixelSpacing(2);
 % % % imgInfo(stackNo).BoundingBox.zMax = double(imgInfo(stackNo).BoundingBox.zMin);
 % % % img{stackNo} = tmp;
 % % % infoPre = info;
 % % %
 % % % numSlicesInStack = [];
 % % % for ii = 2:size(filesList, 1)
 % % %     try
 % % %         info = dicominfo([dicomDir filesList(ii).name]);
 % % %     catch
 % % %         continue;
 % % %     end
 % % %     if ( ~isequal(info.Columns, infoPre.Columns) || ~isequal(info.Rows, infoPre.Rows) || ~isequal(info.PatientOrientation, infoPre.PatientOrientation))
 % % %         img{stackNo} = eval([class(tmp) '(zeros(imgInfo(stackNo).xSize, imgInfo(stackNo).ySize, imgInfo(stackNo).zSize))']);
 % % %         numSlicesInStack(stackNo) = imgInfo(stackNo).zSize;
 % % %         stackNo = stackNo + 1;
 % % %         try
 % % %             imgInfo(stackNo).BoundingBox.xMin = info.ImagePositionPatient(1);
 % % %         catch
 % % %             imgInfo(stackNo).BoundingBox.xMin = 0;
 % % %         end
 % % %         try
 % % %             imgInfo(stackNo).BoundingBox.yMin = info.ImagePositionPatient(2);
 % % %         catch
 % % %             imgInfo(stackNo).BoundingBox.yMin = 0;
 % % %         end
 % % %         try
 % % %             imgInfo(stackNo).BoundingBox.zMin = info.ImagePositionPatient(3);
 % % %         catch
 % % %             imgInfo(stackNo).BoundingBox.zMin = 0;
 % % %         end
 % % %         try
 % % %             imgInfo(stackNo).PixelSpacing     = info.PixelSpacing;
 % % %         catch
 % % %             imgInfo(stackNo).PixelSpacing     = 1;
 % % %         end
 % % %
 % % %         imgInfo(stackNo).SliceThickness   = 0;
 % % %         try
 % % %             imgInfo(stackNo).xSize = info.Columns;
 % % %             imgInfo(stackNo).ySize = info.Rows;
 % % %         catch
 % % %             tmp = dicomread(info);
 % % %             imgInfo(stackNo).xSize = size(tmp, 2);
 % % %             imgInfo(stackNo).ySize = size(tmp, 1);
 % % %         end
 % % %
 % % %
 % % %         imgInfo(stackNo).zSize = 1;
 % % %         imgInfo(stackNo).BoundingBox.xMax = imgInfo(stackNo).BoundingBox.xMin + (imgInfo(stackNo).xSize -1) * imgInfo(stackNo).PixelSpacing(1);
 % % %         imgInfo(stackNo).BoundingBox.yMax = imgInfo(stackNo).BoundingBox.yMin + (imgInfo(stackNo).ySize -1) * imgInfo(stackNo).PixelSpacing(2);
 % % %         imgInfo(stackNo).BoundingBox.zMax = imgInfo(stackNo).BoundingBox.zMin;
 % % %         infoPre = info;
 % % %     else
 % % %         try
 % % %             imgInfo(stackNo).BoundingBox.zMax = info.ImagePositionPatient(3);
 % % %         catch
 % % %             imgInfo(stackNo).BoundingBox.zMax     = imgInfo(stackNo).BoundingBox.zMax + 1;
 % % %         end
 % % %         imgInfo(stackNo).zSize = imgInfo(stackNo).zSize + 1;
 % % %         imgInfo(stackNo).SliceThickness = (imgInfo(stackNo).BoundingBox.zMax - imgInfo(stackNo).BoundingBox.zMin)/imgInfo(stackNo).zSize;
 % % %     end
 % % % end
 % % % img{stackNo} = eval([class(tmp) '(zeros(imgInfo(stackNo).xSize, imgInfo(stackNo).ySize, imgInfo(stackNo).zSize))']);
 % % % numSlicesInStack(stackNo) = imgInfo(stackNo).zSize;
 % % %
 % % % switch opt
 % % %     case 'all'
 % % %         f = 1;
 % % %         for ii = 1:length(numSlicesInStack)
 % % %             for jj = 1:numSlicesInStack(ii)
 % % %                 info = dicominfo([dicomDir filesList(f).name]);
 % % %                 tmp  = dicomread(info);
 % % %                 img{ii}(:, :, jj) = tmp;
 % % %                 f = f + 1;
 % % %             end
 % % %         end
 % % %     case 'first'
 % % %         for ii = 1:numSlicesInStack(1)
 % % %             info = dicominfo([dicomDir filesList(ii).name]);
 % % %             tmp  = dicomread(info);
 % % %             img{1}(:, :, ii) = tmp;
 % % %         end
 % % %         imgInfo = imgInfo(1);
 % % %         img = img{1};
 % % %     case 'last'
 % % %         if numel(numSlicesInStack) == 1
 % % %             fileInd = 0;
 % % %         else
 % % %             fileInd = sum(numSlicesInStack(end-1));
 % % %         end
 % % %         for ii = 1:numSlicesInStack(end)
 % % %             info = dicominfo([dicomDir filesList(ii+fileInd).name]);
 % % %             tmp  = dicomread(info);
 % % %             img{end}(:, :, ii) = tmp;
 % % %         end
 % % %         imgInfo = imgInfo(end);
 % % %         img = img{end};
 % % %     otherwise
 % % %         [~, ind] = max(numSlicesInStack);
 % % %         if ind == 1
 % % %             fileInd = 0;
 % % %         else
 % % %             fileInd = sum(numSlicesInStack(1:ind-1));
 % % %         end
 % % %         for ii = 1:numSlicesInStack(ind)
 % % %             info = dicominfo([dicomDir filesList(ii+fileInd).name]);
 % % %             tmp  = dicomread(info);
 % % %             img{ind}(:, :, ii) = tmp;
 % % %         end
 % % %         imgInfo = imgInfo(ind);
 % % %         img = img{ind};
 % % % end