How to convert stance image from index to grayscale

I am looking to replace the 'ind2gray' function as we no longer have access to the image processing toolbox. Is there a script to replace it?

% Convert stance image for each phase from 'index' to 'grayscale' and
% insert into Phases field: ManipulatedStances
for k = 1 : numPhases;
    thisPhase = (Phases(k).Stance);
    I = ind2gray(thisPhase,map1);
    Phases(k).ManipulatedStances = I; 
end

Answer 1

[DGM]

No, don't use BT709 luma.

No, don't discard columns from the colormap.

If you're using IPT ind2gray(), just use ind2rgb(). It's not part of IPT, nor was it in IPT in 2016.

% get an indexed color image
[indpict CT] = imread('canoe.tif');

% convert to RGB
rgbpict = ind2rgb(indpict,CT);

% calculate luma
% since we don't have IPT, and we're living in early 2016, do it manually
factors = permute([0.299 0.587 0.114],[1 3 2]); % Rec 470/601 (this is what ind2gray would use)
ypict = sum(bsxfun(@times,im2double(rgbpict),factors),3);
ypict = im2uint8(ypict); % assuming you want uint8 output

% is there a difference? no, not really.
ypict0 = ind2gray(indpict,CT); % for comparison
immse(ypict0,ypict) % nothing but rounding errors (1LSB in a few specks)

ans = 0.0017

imshow(ypict)

Answer 2

[John BG]

Hi Denni

You want to go from RGB to Y. Y is called luminance.

There is a good presentation of TV colour calibration charts in

http://bennettcain.com/blog/2012/2/21/hd-monitor-calibration-smpte-arib-and-white-balance.html

Basically, there are 2 industry accepted standards to calculate Y from RGB:

A=imread('PAL TV colour calibration chart 1.jpg')
imshow(A)

1.-

A1=A(:,:,1);A2=A(:,:,2);A3=A(:,:,3)
Y_CCIR601=.299*A1+.587*A2+.114*A3
imwrite(Y_CCIR601,'Y_CCIR601.png')  % save BW image in graphics file format

2.-

Y_ITUR=.2126*A1+.7152*A2+.0722*A3
imwrite(Y_ITUR,'Y_ITUR.png')  % save BW image in graphics file format

If you find this answer of any help solving your question, please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]

Answer 3

[Image Analyst]

If you have a map, that converts indexed into gray, you can do it with a for loop. If the map is a 256 by 3, then extract just the first column so it's a 1-D vector of 256 elements.

if size(map, 2) > 1
    map = map(:, 1);
end

Since you don't have intlut() either (because you no longer have the toolbox), you'll have to use a for loop

grayImage = zeros(size(thisPhase)); % Initialize
for k = 1 : numel(grayImage)
    inputIndex = thisPhase(k);
    grayImage(k) = uint8(255 * map(inputIndex));
end

Now grayImage will be a uint8 2D array the same size as thisPhase with gray levels in the range 0-255 that were assigned according to the map that you appear to have already.

[Denni Purcell]

I have altered the script to look like this, however I am receiving an error code

for k = 1 : numPhases; thisPhase = (Phases(k).Stance); grayImage = zeros(size(thisPhase)); % Initialize end

for k = 1 : numel(grayImage) inputIndex = thisPhase(k); grayImage(k) = uint8(map(inputIndex)); end

Error code Undefined function or variable 'map'.

Error in Matlab_PBG_code (line 208) grayImage(k) = uint8(map(inputIndex));

[Image Analyst]

Sorry - my typo. According to your original code, your variable was called map1 instead of map, so change it to map1 and it should work.

[Denni Purcell]

Sorry still having errors...

Undefined function or variable 'map1'.

Error in Matlab_PBG_code (line 209)
    grayImage(k) = uint8(map1(inputIndex));

Previously in the code this is how the map1 is defined

map = colormap(jet(pCap)); map(1,:) = [1 1 1];
colormap(map);

Any suggestions here?

[Image Analyst]

This was your code:

I = ind2gray(thisPhase,map1);

so where did you get map1???? I have no idea because you did not show me that code. It would have to be something or else you couldn't even do what you showed, even if you did have the toolbox. So show the code from that program where map1 was first assigned.

[Denni Purcell]

sorry, this code was not written by me and I have minimal knowledge of Matlab. This is the code from the start where the map is labeled

%%--------------------------------------------------------------------- %%
% ------------------------- 01: LOAD DATA FILES ------------------------- %
% ----------------------------------------------------------------------- %
tic 

clear all; close all; clc

% Patient Name
% prompt1 = {'First Name, Last Name -- DD/MM/YY (Date of Assessment)'};
% dlg_title1 = 'Patient Assessment Details';
% answer1 = inputdlg(prompt1,dlg_title1, [1, 65]);
%                                                                           
% PatInf = answer1;

% Enter session path
prompt = {'Enter patient 3DGA folder path'};
dlg_title = 'Input path';
pname = inputdlg(prompt, dlg_title);
pname = [char(pname) ''];

cd(pname)

% COF File
COF_File = uigetfile('.csv','Select Centre Of Force Data File');
    % Footstrike Data File
Footstrike_File = uigetfile('.csv','Select Force Data All File');
    % Frame Data File
Frame_Data_File = uigetfile('.csv','Select Frame Data File');
    % Stance File
Stance_File = uigetfile('.csv','Select Stance Data File');
                                                                           % IMPORTANT!!!
    % "Pressure cap" (for the pressure-colour scale)
prompt = {'Enter the pressure threshold'};
name = 'Pressure cap';
defaultans = {'50'};
answer = inputdlg(prompt,name,[1 20],defaultans);
pCap = str2num(answer{:});

% Set custom colormap
map = colormap(jet(pCap)); map(1,:) = [1 1 1];
colormap(map);
% IMPORTANT!!!
% If the figure's show that the pCap has been set too low or too high,
% leave the figures open and change the value for pCap to whatever. Run
% lines 28 to 31 to update the pressure map (you can do this by
% highlighting those lines then pressing 'F9'). If only one of the figure
% windows updates, then you'll need to repeat this for each figure window.
% Simply open an unaffected figure window, then switch back to the editor
% tab (this one) and run lines 28 to 31.
% IMPORTANT!!!

% Defining the colormap without a current figure window means MATLAB will
% generate an empty one. The folowing line closes it.
close all;

      % Determine the number of frames
  nFrames = length(xlsread(COF_File));

Then this section is where the I = ind2gray(thisPhase,map1); is used

   %%--------------------------------------------------------------------- %%
  % -------------------------- 02: COMPILE DATA --------------------------- %
  % ----------------------------------------------------------------------- %

% ---- STRUCTURE: Phases ------------------------------------------------ %
% This structure contains most of the information of the loaded data: the
% number of phases, how many strikes per phase, the cropping information
% for each footstrike and more... check it out once this section is run!

% Read entire stance spreadsheet from the first numerical cell to the last
Stance = xlsread(Stance_File);

% Find number of rows (and cols) that were extracted into "Stance"
[nRows,~]=size(Stance);

exRows = nRows + 42;

% Now extract only the data of interest from spreadsheet
% NEW_Stance = xlsread(Stance_File,strcat('A42',:,'CJ',num2str(exRows))); %
NEW_Stance = xlsread(Stance_File, ['A42:CJ' num2str(exRows)]);

% Determine the number of phases in the files selected
% ('+2' because the top phase never has the top two rows to begin with, 
% unlike subsequent rows. '/386' because there are 386 rows in each phase 
% (resolution plus two non-numerical rows - open the .csv to see what this 
% means))
numPhases = (nRows + 2)/386;

% Create blank two rows of non-numbers, and add them to the data extracted
% into NEW_Stance
two_rows = ones(2,88);
two_rows(two_rows==1)=NaN;
NEW_Stance = [two_rows;NEW_Stance];

% Create empty structure with many fields.
Phases = struct('Stance',{},'COF',{},'ManipulatedStances',{},'Feet_Props',{},'Sorted_BB',{},'NumberOfStrikes',{},'Feet_Crops',{});

% Separate each phase from each other
for k = 1 : numPhases;
    % r1 and r2 determine the first and last row for the specific phase
    r1 = ((k-1)*386)+3;
    r2 = ((k-1)*386)+386;
    % here, NEW_Stance is cropped from r1 to r2, and columns 1 to 88. This
    % crop defines one stance, and it is placed into its own cell in the
    % Stance field of the structure Phases
    Phases(k).Stance = NEW_Stance(r1:r2,1:88);
end

% Read the COF file
% The numbers ",38,3,[38 3 (37+nFrames) 4])" define where MATLAB will
% extract from, out of the COF file. Open the .csv file to see that the
% first 38 lines is gibberish.
COF = csvread(COF_File,38,3,[38 3 (37+nFrames) 4]); 

% Convert the COF data into binary, and single out only one column. This
% makes segmentation into phases easier
bwCOF=COF;
bwCOF(bwCOF==-1)=0;
bwCOF(bwCOF>0)=1;
bwCOF=bwCOF(:,1);

% Segmentation into phases
bwCOF=[0;bwCOF;0];
D=diff(bwCOF);
b.start = find(D == 1);
b.finish = find (D == -1) -1;

% Compilation of the phases into a single matrix
b.start=transpose(b.start); 
b.finish=transpose(b.finish);
Dd = transpose([b.start; b.finish]);
% You have with this matrix, Dd, the start and end rows, for each phase!!!
                                                                           % IMPORTANT!!!
% IMPORTANT!!!
% Opening the COF .csv file shows that after some number of rows, there are
% a small bunch of rows with '-1' in them, after which the real data
% continues. This indicates where one phase ends and another begins, and
% this is how the foregoing lines determined where each phase began and
% ended. This is also why it is important to cut the frames exactly as the
% phases are defined in Walkway (as opposed to adding 6 and taking 6
% frames); doing so will confuse the COF data. Check the data and make sure
% that there isn't just a few lines of real data followed by the -1s. The 
% data should either start with the -1s, or start with many rows of real 
% data. If there is only one or a few lines of real data at the beginning,
% then the frames should either be re-cut from Walkway to make sure it was
% done correctly, OR the first few lines should be edited and changed to
% -1.
% IMPORTANT!!!

% Use a for-loop to fill the 'COF' field with each separate phase's COF data
for k = 1 : numPhases;
    start = Dd(k,1);
    finish = Dd(k,2); 
    Phases(k).COF = COF(start:finish,1:2);
end

% Convert stance image for each phase from 'index' to 'grayscale' and
% insert into Phases field: ManipulatedStances
for k = 1 : numPhases;
    thisPhase = (Phases(k).Stance);
    I = ind2gray(thisPhase,map1);
    Phases(k).ManipulatedStances = I; 
end

% Convert the now grayscale images to binary.
thresholdValue = 1;
for k = 1 : numPhases;
    thisPhase = (Phases(k).ManipulatedStances);
    I = thisPhase < thresholdValue;
    Phases(k).ManipulatedStances = I; 
end
                                                                           % IMPORTANT!!!
% Now patch the images so each foot strike becomes a single and separate 
% region
se = strel('square',25);
for k = 1 : numPhases;
    thisPhase = (Phases(k).ManipulatedStances);
    I = imclose(thisPhase,se);
    Phases(k).ManipulatedStances = I; 
end
% IMPORTANT!!!
% Before lines 166-173, what MATLAB sees is exactly what we see, which
% means that there is some space between the toes and the feet, and the
% heel is often separate as well. MATLAB would then think that each toe and
% heel are separate footstrikes. Lines 166-173 attempt to join them
% together, so that they are contiguous (no zeros between) elements. It
% does this with a patching element, defined by the variable 'se'. By
% default this is a square shaped element with sides 25 elements long. This
% section can be edited to change the size of the square, or the shape of
% the element as well if the feet are not cropped out correctly. Type 'help
% strel' or 'doc strel' to see more information on patching elements.
% IMPORTANT!!!

% Determine the bounding boxes of each phase's footstrikes, and how many 
% there are in each phase
for k = 1 : numPhases;
    thisPhase = (Phases(k).Stance);
    thisManipulatedPhase = (Phases(k).ManipulatedStances);
    FeetProps = regionprops(thisManipulatedPhase, thisPhase,'BoundingBox');
    Phases(k).Feet_Props = FeetProps;
    Phases(k).NumberOfStrikes = size(Phases(k).Feet_Props, 1);
end
% NOTE: 
% When displayed (in the rep. feet section), the number of rows will
% correspond to the number of phases, and the number of columns will
% correspond to the number of footstrikes in the phase which has the MOST
% (i.e. if three out of four phases have three strikes, but one of them has
% five, there will be five columns).

% Order the bounding boxes to correspond with the actual footsrikes
for k1 = 1 : numPhases;
    for k = 1 : Phases(k1).NumberOfStrikes;
        thisStrikesBoundingBox = Phases(k1).Feet_Props(k).BoundingBox;
        Phases(k1).Sorted_BB(k,:) = thisStrikesBoundingBox;
    end
    clear k;
end

for k1 = 1 : numPhases;
    Phases(k1).Sorted_BB = sortrows(Phases(k1).Sorted_BB,2);
    thisPhasesCOF_start = Phases(k1).COF(1);
    thisPhasesCOF_end   = Phases(k1).COF(end,1);
        if thisPhasesCOF_start < thisPhasesCOF_end;
        else
            Phases(k1).Sorted_BB = flip(Phases(k1).Sorted_BB);
        end
    for k = 1 : Phases(k1).NumberOfStrikes;
        BoundingBox = Phases(k1).Sorted_BB(k,:);
        Phases(k1).Feet_Props(k).OrderedBoundingBoxes = BoundingBox;
    end
end

% Input data of each cropped foot strike into the Phases field,
% "Feet_Crops"

% Determine the number of strikes in the phase which has the most of them.
for k = 1 : numPhases
    maxStrikes(k,:) = Phases(k).NumberOfStrikes;
end
MaxNoStrikes = max(maxStrikes);

for k1 = 1 : numPhases;
    for k = 1 : Phases(k1).NumberOfStrikes;
        thisPhase = Phases(k1).Stance;
        thisStrikesBoundingBox = Phases(k1).Feet_Props(k).OrderedBoundingBoxes;
        Phases(k1).Feet_Crops(k).CroppedImages = imcrop(thisPhase,thisStrikesBoundingBox);
        kk = num2str(k+MaxNoStrikes*(k1-1));
        Phases(k1).Feet_Crops(k).Name = ['strike no.', char(kk)];
    end
end

% Flip footstrikes so that they are all the same orientation. They will now
% either all be upside-down, or all the right way up.
for k1 = 1 : 2 : numPhases; % Odd numbered phases
    for k = 1 : Phases(k1).NumberOfStrikes;
        Phases(k1).Feet_Crops(k).CroppedImages = flip(Phases(k1).Feet_Crops(k).CroppedImages);
        Phases(k1).Feet_Crops(k).CroppedImages = fliplr(Phases(k1).Feet_Crops(k).CroppedImages);
    end
end

% Create vector listing all the strikes
str = {}; 
% str = cell(1,< TOTAL NUMBER OF STRIKES >); 
% ^^ Use once you determine how to call the total number of strikes
for k = 1 : numPhases;
    str{k} = {Phases(k).Feet_Crops.Name};
end
str = horzcat(str{:});

% ---- STRUCTURE: All_Feet_Crops ---------------------------------------- %
% This is a simple structure with only one field. It each cell holds the
% image of each footstrike. Used in section 04.

% Concatenate all footstrikes' bounding boxes into a single, cell array
All_Feet_Crops = {}; 
% All_Feet_Crops = cell(1,< TOTAL NUMBER OF STRIKES >);
% ^^ Use once you determine how to call the total number of strikes
for k = 1 : numPhases;
    All_Feet_Crops{k} = {Phases(k).Feet_Crops.CroppedImages};
end
All_Feet_Crops = horzcat(All_Feet_Crops{:});

% ---- STRUCTURE: All_Feet_Crops_GOOD ----------------------------------- %
% This structure has more information on each footstrike than
% All_Feet_Crops, including the frame range for each strike. Used in
% section 05.

% Read Foot Strike File
Footstrike_Data = xlsread(Footstrike_File);
% Determine size
[nRows,nCols]=size(Footstrike_Data);
% Now extract only the relevant information
Footstrike_Data = csvread(Footstrike_File,38,3,[38 3 (37+nFrames) (nCols-1)]);

% Re-order the Footstrike_Data columns to reflect the order of the strikes
[fStrikes_rows,fStrikes_cols] = size(Footstrike_Data);
Footstrike_Data = [zeros(1,fStrikes_cols); Footstrike_Data];
for k = 1 : fStrikes_cols; % remember to change 1 back to 4
    Footstrike_Data(1,k) = find(Footstrike_Data(:,k),1);
end
[YY,II] = sort(Footstrike_Data(1,:));
Footstrike_Data = Footstrike_Data(:,II);
Footstrike_Data(1,:) = [];

% Create structure that lists all the cropped feetstrikes, and their
% corresponding frame ranges

% List of all feet crops
All_Feet_Crops2 = struct('Crops',{},'First_Frame',{},'Last_Frame',{});
for k1 = 1 : numPhases;
    for k = 1 : Phases(k1).NumberOfStrikes;
        ThisCroppedImage = Phases(k1).Feet_Crops(k).CroppedImages;
        All_Feet_Crops2(k+MaxNoStrikes*(k1-1)).Crops = ThisCroppedImage;
    end
end

% Remove any possible empty spaces from feetstrike structure, by moving
% data into a new structure
count = 1;
for f=1:length(All_Feet_Crops2)
    if(~isempty(All_Feet_Crops2(f).Crops))
        All_Feet_Crops_GOOD(count) = All_Feet_Crops2(f);
        count = count + 1;
    end
end

% Delete useless structure "All_Feet_Crops2"
clear All_Feet_Crops2;

% Extract each footstrike's frame range
for k = 1 : fStrikes_cols;
    ThisStrikesFirstFrame = find(Footstrike_Data(:,k),1,'first');
    ThisStrikesLastFrame = find(Footstrike_Data(:,k),1,'last');
    All_Feet_Crops_GOOD(k).First_Frame = ThisStrikesFirstFrame;
    All_Feet_Crops_GOOD(k).Last_Frame = ThisStrikesLastFrame;
end

% Add the Bounding Box information to each footstrike
count = 0;
for k1 = 1 : numPhases;
    for k = 1 : Phases(k1).NumberOfStrikes;
        ThisStrikesBoundingBox = Phases(k1).Feet_Props(k).OrderedBoundingBoxes;
        All_Feet_Crops_GOOD(count + 1).BoundingBox = ThisStrikesBoundingBox;
        count = count + 1;
    end
end

% Add footstrike number to each footstrike
for k = 1 : length(All_Feet_Crops_GOOD);
    kk = num2str(k);
    All_Feet_Crops_GOOD(k).Number = ['strike no.', char(kk)];
end

% Create vector listing all the strikes
str2 = cell(1,length(All_Feet_Crops_GOOD));
for k = 1 : length(All_Feet_Crops_GOOD);
    str2{k} = {All_Feet_Crops_GOOD(k).Number};
end
str2 = horzcat(str2{:});

So now I have replaced the ind2gray section with this

for k = 1 : numPhases;
    thisPhase = (Phases(k).Stance);
    grayImage = zeros(size(thisPhase)); % Initialize
    end

for k = 1 : numel(grayImage)
    inputIndex = thisPhase(k);
    grayImage(k) = uint8(map1(inputIndex));
end

And have been getting the errors mentioned in the previous post

[Image Analyst]

I have no idea where they got map1. It should have thrown an error since it's not defined before they started using it.