function extractnodules( )

%EXTRACTNODULES Extract nodule images from lung CT scans

%

%   Ekarin Varutbangkul

%   July 8, 2006

%

%   Reads nodule data from 'DCMInfo2.txt' in working directory and pulls

%   individual nodule images from DICOM images referenced in the data

%   file, writing them to new DICOM images in the 'nodules/' directory.

%   Also saves all Haralick and annotation data to 'haralick.txt' in the

%   same directory.

tic

fid = fopen('contours.txt', 'rt');

% read line from data file

tline = fgetl(fid);

% extract nodule information

[pk n_filename n_minx n_maxx n_miny n_maxy n_x n_y ] ...

     = strread(tline, '%d %s %d %d %d %d %d %d');

%cn = 1;

%while cn<4 

while feof(fid) == 0

%    k1 = n_file;

%    k2 = n_id;

    k1 = pk;

    s_minx = n_minx;

    s_maxx = n_maxx;

    s_miny = n_miny;

    s_maxy = n_maxy;

    fn = sprintf('%s', char(n_filename));

    fn2 = sprintf('%d.dcm', pk);

    fprintf(1,'    Processing: %d  (x=[%d,%d] y=[%d,%d]) ...\n', pk, n_minx, n_maxx, n_miny, n_maxy);

    % make sure file exists

    if (exist(fn,'file'))

        % read image

        info = dicominfo(fn);

        I = dicomread(info);

        % get nodule range

        n_rows = n_maxy - n_miny + 1;

        n_cols = n_maxx - n_minx + 1;

        % minimum nodule size is 5x5 pixels

        if (n_rows > 4 && n_cols > 4)

            % extract nodule data

            I2 = int16(zeros(n_rows,n_cols));

            I3 = int16(zeros(n_rows,n_cols));

            I4 = int16(zeros(n_rows,n_cols));

            for r = 1:n_rows

                for c = 1:n_cols

                    I2(r,c) = I(n_miny+(r-1),n_minx+(c-1));

                end

            end                

            bg = -2000;

            k_minx = n_minx;

            k_miny = n_miny;

           while (k1 == pk) & (feof(fid) == 0)

                I3((n_y-n_miny)+1,(n_x-n_minx)+1)= 1;  

                if feof(fid) == 0

                    tline = fgetl(fid);

                    [pk n_filename n_minx n_maxx n_miny n_maxy n_x n_y ] ...

                    = strread(tline, '%d %s %d %d %d %d %d %d');   

                end

           end

           BW = zeros(n_rows,n_cols);

           BW = imfill(I3,'holes');

           I3 = BW;

           I3(I3 == 0) = bg;

           I4 = BW;

% old segmentation code 

%{ 

            for r = 1:n_rows

                for c = 1:n_cols

                    if (I3(r,c) ~= 1)

                        I3(r,c) = bg;

                    else

                        break;

                    end

                end

            end

            for r = 1:n_rows

                for c = n_cols:-1:1

                    if (I3(r,c) ~= 1)

                        I3(r,c) = bg;

                    else

                        break;

                    end

                end

            end

            for c = 1:n_cols

                for r = 1:n_rows

                    if (I3(r,c) ~= 1)

                        I3(r,c) = bg;

                    else

                        break;

                    end

                end

            end

            for c = 1:n_cols

                for r = n_rows:-1:1

                    if (I3(r,c) ~= 1)

                        I3(r,c) = bg;

                    else

                        break;

                    end

                end

            end

%}

            for r = 1:n_rows

                for c = 1:n_cols

                    if I3(r,c)~= bg

                        I3(r,c) = I(k_miny+(r-1),k_minx+(c-1));

                    end

                end

            end 

            for r = 1:n_rows

                for c = 1:n_cols

                    if I4(r,c)== 0

                        I4(r,c) = I(k_miny+(r-1),k_minx+(c-1));

                    else

                        I4(r,c) = bg;

                    end

                end

            end 

            % save nodule image

            info.Width = n_cols;

            info.Height = n_rows;

            info.Rows = n_rows;

            info.Columns = n_cols;

            dicomwrite(I2, strcat('crop/',fn2),info);

            dicomwrite(I3, strcat('segmented_nodules/',fn2),info);

            dicomwrite(I4, strcat('segmented_background/',fn2),info);

%            fprintf(1,'    Saving data ...\n');

         else

%            fprintf(1,'%s is too small (less than 5x5) -- skipping ...\n',fn2)

            epk = pk;

            while (feof(fid) == 0) & (epk == pk)        

                tline = fgetl(fid);

                    [pk n_filename n_minx n_maxx n_miny n_maxy n_x n_y ] ...

                    = strread(tline, '%d %s %d %d %d %d %d %d');   

            end

        end

    else

        fprintf(1,'    File does not exist: %s -- skipping ...\n', fn)

        while (feof(fid) == 0) & (epk == pk)  

            tline = fgetl(fid);

                    [pk n_filename n_minx n_maxx n_miny n_maxy n_x n_y ] ...

                    = strread(tline, '%d %s %d %d %d %d %d %d');   

        end

    end

%    cn = cn+1;

end

fclose(fid);

toc