Here are the programs that were developped to achieve the automatic computation of the cell size. You probably don't need to be very familiar with Matlab to read them and understand how they work, and the comments (announced by the sign %) are here to help!

The following links will bring you directly to the program codes:

cellsizer.m......... main program

filtr.m................ filtering program

matlabsobel.m... Matlab-based edge enhancement program

persobel.m........ home-made edge enhancement program

calcpsd.m......... PSD computation

disppsd.m......... PSD display

angplot.m.......... angular scanning

radplot.m.......... radial scanning

peak.m.............. ultimate cell size computation

center.m............ detection of the center of the PSD

grayscale.m....... picture grayscaling

Main program: cellsizer.m(download)

function RESULT=cellsizer(MATRIX)

% Enters a serie of procedures to process

% the cellular pattern analysis.

% You can display the cellular pattern, visualize the PSD,

% launch an averaged angular scanning, which will help you

% defining the angle of the PSD legs, and use it in an

% averaged radial scanning to extract the spatial wave length

% of the cellular pattern.

% This programs returns an array:

% 1st column: cell size harmonics

% 2nd column: relative energy of each harmonic

% The picture must be Black & White.

% IT IS ALSO STRONGLY RECOMMENDED, IF NOT REQUIRED, TO USE SQUARE IMAGES.

% NOT DOING SO CAN LEAD TO FALSE RESULTS...

% If possible, use an image whose dimensions in pixels are

% a power of 2.

% EXAMPLES:

% >> cellsizer('cell.jpg');

% where cell.jpg is the cellular pattern image

% OR

% >> cellsizer(MATRIX);

% where MATRIX is the double-type matrix of

% the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg')).

% To keep the harmonics & associated energies in the variable ARRAY, type

% >> ARRAY=cellsizer(...);

if isstr(MATRIX)

MATRIX=double(imread(MATRIX)); % Turns the image into a matrix

end

pi=3.14159265;

[HEIGHT WIDTH]=size(MATRIX);

if HEIGHT~=WIDTH

disp(' ')

disp(' WARNING: Your image isn''t square.');

disp(' You may love taking risks, or just getting wrong results...');

end

FINAL=calcpsd(MATRIX,HEIGHT,WIDTH); % DSP calculation

[Xc Yc]=center(FINAL);

m=0; % Default value for the main menu

while m~=7 % If m==7 : QUIT

m=menu('WHAT WOULD YOU LIKE TO DO?','DISPLAY CELLULAR PATTERN AND PSD REPRESENTATION','FILTERING','EDGE ENHANCEMENT','ANGULAR SCANNING','RADIAL SCANNING','CELL WIDTH','QUIT and learn why I wrote this program...');

if m==1

disppsd(MATRIX,FINAL); % DSP visualization

end

if m==2

FILTERED=filtr(MATRIX,HEIGHT,WIDTH,Xc,Yc); % Image filtering

sm=menu('Would you like to work on the filtered image as of now?','Yes','No');

if sm==1

MATRIX=double(FILTERED);

FINAL=calcpsd(MATRIX,HEIGHT,WIDTH); % DSP calculation for the filtered image

end

if m==3 % Edge enhancement

sm=menu('Which program do you wish to use?','Matlab program','Home-made program');

if sm==1, EDGED=matlabsobel(MATRIX);

else, EDGED=persobel(MATRIX); end

ssm=menu('Would you like to work on the edge enhanced image as of now?','Yes','No');

if ssm==1

MATRIX=double(EDGED);

FINAL=calcpsd(MATRIX,HEIGHT,WIDTH); % DSP calculation for the filtered image

end

if m==4

angplot(FINAL,HEIGHT,WIDTH,Xc,Yc); % Angular scanning

end

if m==5

ANGL=radplot(FINAL,HEIGHT,WIDTH,Xc,Yc); % Radial scanning. Returns the angle of the scanned leg, to use it in peak.

end

if m==6

if HEIGHT==WIDTH

RESULT=peak(HEIGHT,ANGL); % Cell width calculation

else

disp(' Impossible to measure the cell size of a non-square image')

end

disp(' ')

disp(' Thanks for running me!')

disp(' Hope to see you again soon...')

disp(' ')

disp('In fact, the reason why I wrote this program is')

why

disp(' ')

Filtering program: filtr.m(download)

function FILTERED=filtr(MATRIX,HEIGHT,WIDTH,Xc,Yc)

% Erases selected areas of the PSD.

% Considering the PSD visualization and its 2 two caracteristics legs,

% you can define 4 areas : 3, 6, 9, 12 o'clock.

% filtr allows you to erase 1 or 2 pairs of opposite areas,

% staying remote of a distance gap/2 (that you are invited to input)

% from the legs you want to keep intact.

% EXAMPLES:

% >> filtr('cell.jpg')

% where cell.jpg is the cellular pattern image.

% OR

% >> filtr(MATRIX)

% where MATRIX is the double-type matrix of

% the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg'));

% You can also specify the HEIGHT, the WIDTH, and the center coordinates Xc and Yc

% of your image. If you do so, type:

% >> filtr(...,HEIGHT,WIDTH,Xc,Yc)

% Enter either 1 argument (the image or its matrix), or 5.

% >> FILTERED=filtr(...) returns as FILTERED the filtered cellular pattern matrix

% (CAUTION: it's not the PSD matrix!).

if isstr(MATRIX)

MATRIX=double(imread(MATRIX)); % Turns the image into a matrix

end

% The previous loop allows us to input the image either as a function argument

% (then CELLIMAGE is a string, for instance cell.jpg),

% or the image matrix (then CELLIMAGE is a matrix).

if nargin==1

FINAL=calcpsd(MATRIX);

[HEIGHT WIDTH]=size(FINAL);

[Xc Yc]=center(FINAL);

end

if nargin~=1 & nargin~=5

error (' Number of expected input arguments: 1 or 5');

end

% The difficulty is to paint black at the right location of the PSD,

% but we can't do it on the actual PSD, since it is then impossible

% to come back to the image matrix (because of .^2, log,... operators).

% So we process the filter on the NONFINAL matrix, whose pixels are

% on the same location as in the PSD (with of course not the same intensity),

% but which allows us to come back to the image afterwards.

m=0;

WIND=MATRIX.*fspecial('gaussian',[HEIGHT WIDTH],0.2*min(HEIGHT,WIDTH)); % Windows the image with a gaussian filter

WIND=grayscale(WIND); % Rescales the entries between 0 and 255.

NONFINAL=fftshift(fft2(WIND));

while m~=2

disp(' ');

disp(' Please enter the angles of the 2 legs of the PSD :');

TH1=2*pi/360*input('\n Positive angle (in degrees) = ');

TH2=2*pi/360*input('\n Negative angle (in degrees) = ');

g=input('\n Surviving gap = ');

g1=g/2/cos(TH1);

g2=g/2/cos(TH2);

sm=menu('What area would you like to erase?','Around 3 and 9 o''clock', 'Around 6 and 12 o''clock','Around 3,6,9 and 12 o''clock');

MASK=ones(HEIGHT,WIDTH);

I=0;

h=waitbar(0,'Please wait, filtering in progress...');

for Y=1:WIDTH

waitbar(I/WIDTH,h)

I=I+1;

X1=Xc-(Y-Yc)*tan(TH1);

X2=Xc-(Y-Yc)*tan(TH2);

for X=1:HEIGHT

if (sm==1 | sm==3) & ((Y<=Yc & X2+g2<X & X<X1-g1) | (Y>=Yc & X1+g1<X & X<X2-g2)), MASK(X,Y)=0; end

if (sm==2 | sm==3) & ((X<X2-g2 & X<X1-g1) | (X1+g1<X & X2+g2<X)), MASK(X,Y)=0; end

end

close(h)

MASKED=MASK.*NONFINAL;

FILTERED=real(ifft2(ifftshift(MASKED)));

LOCALFINAL=grayscale(calcpsd(FILTERED));

figure;colormap(gray);imshow(LOCALFINAL);image(LOCALFINAL);title('Filtr -- Filtered Power Spectral Density');axis on;zoom on

figure;colormap(gray);imshow(FILTERED);image(FILTERED);title('Filtr -- Filtered Cellular Pattern');axis on;zoom on % Imshow>>Image gives better results than Image alone. I don't know why...

m=menu('What would you like to do now?','Try again','Quit or return to main menu');

end

return

Matlab-based edge enhancement program: matlabsobel.m(download)

function EDGED=matlabsobel(MATRIX)

% Edge enhancement using the Sobel method.

% EXAMPLES:

% >> matlabsobel('cell.jpg')

% where cell.jpg is the cellular pattern image

% OR

% >> matlabsobel(MATRIX)

% where MATRIX is the double-type matrix of

% the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg')).

% >> EDGED=matlabsobel(...) returns as EDGED the matrix of the

% enhanced edge image.

if isstr(MATRIX)

MATRIX=double(imread(MATRIX)); % Turns the image into a matrix

end

% The previous loop allows us to input the image either as a function argument

% (then CELLIMAGE is a string, for instance cell.jpg),

% or the image matrix (then CELLIMAGE is a matrix).

[EDGED,t]=edge(MATRIX,'sobel'); % First try: Matlab selects its own threshold.

EDGED=grayscale(double(EDGED));

figure;colormap(gray);imshow(EDGED);title(sprintf('Matlabsobel -- Enhanced edges with a threshold of: %d',t));axis on;

m=menu('What would you like to do now?','Change threshold','Quit or return to main menu');

while m~=2

sprintf(' Previous used threshold: %d',t)

t=input(' Enter new threshold: ');

EDGED=grayscale(double(edge(MATRIX,'sobel',t)));

figure;colormap(gray);imshow(EDGED);title(sprintf('Matlabsobel -- Enhanced edges with a threshold of: %d',t));axis on;

m=menu('What would you like to do now?','Change threshold','Quit or return to main menu');

end

return

Home-made edge enhancement program: persobel.m(download)

function EDGED=persobel(MATRIX)

% Personnal version of edge enhancement,

% using the Sobel method.

% EXAMPLES:

% >> persobel('cell.jpg')

% where cell.jpg is the cellular pattern image

% OR

% >> persobel(MATRIX)

% where MATRIX is the double-type matrix of

% the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg')).

% >> EDGED=persobel(...) returns as EDGED the matrix of the

% enhanced edge image.

if isstr(MATRIX)

MATRIX=double(imread(MATRIX)); % Turns the image into a matrix

end

%Mv=[-1 0 1 ; -2 0 2 ; -1 0 1];

%Mh=[1 2 1 ; 0 0 0 ; -1 -2 -1];

[HEIGHT WIDTH]=size(MATRIX);

GRAD=zeros(HEIGHT,WIDTH);

EDGED=zeros(HEIGHT,WIDTH);

%figure;imshow(MATRIX);colormap(gray);title('Cellular Pattern');axis on;

%zoom on

h=waitbar(0,'Computing the image gradient, patience required!');

for I=2:HEIGHT-1 % This loop computes the gradient matrix GRAD.

waitbar(I/HEIGHT,h)

for J=2:WIDTH-1

GRAD(I,J)=(MATRIX(I-1,J+1)+2*MATRIX(I,J+1)+MATRIX(I+1,J+1)-MATRIX(I-1,J-1)-2*MATRIX(I,J-1)-MATRIX(I+1,J-1))^2+(MATRIX(I-1,J-1)+2*MATRIX(I-1,J)+MATRIX(I-1,J+1)-MATRIX(I+1,J-1)-2*MATRIX(I+1,J)-MATRIX(I+1,J+1))^2;

end

close(h)

M=max(max(GRAD));

m=1;

while m~=2 % This loop compares the gradient to the input threshold.

t=input('\n Enter threshold percentage (for instance 20): ');

h=waitbar(0,'Comparing the gradient to the threshold, so far so good...');

for I=2:HEIGHT-1

waitbar(I/HEIGHT,h)

for J=2:WIDTH-1

if GRAD(I,J)>t*M/100

EDGED(I,J)=255;

end

close(h)

figure;colormap(gray);image(EDGED);

title(sprintf('Persobel -- Edge enhancement -- Threshold : %d %%',t));

axis on;zoom on;

m=menu('What would you like to do now?','Change threshold','Quit or return to main menu');

end

PSD computation: calcpsd.m(download)

function FINAL=calcpsd(MATRIX,HEIGHT,WIDTH);

% Computes the DSP calculation and returns the matrix of the DSP.

% EXAMPLES:

% >> calcpsd('cell.jpg')

% where cell.jpg is the cellular pattern image

% OR

% >> calcpsd(MATRIX)

% where MATRIX is the double-type matrix of

% the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg')).

% You can also specify the dimensions of the image/matrix

% (they are the same), but it's not compulsory.

% In this case, type:

% >> calcpsd(...,HEIGHT,WIDTH)

% where HEIGHT and WIDTH are the number of rows and columns.

% If you want to get the PSD matrix and call it, say, FINAL, type:

% >> FINAL=calcpsd(...);

% Don't forget the semi-colon, or the whole matrix will be displayed!

if isstr(MATRIX)

MATRIX=double(imread(MATRIX)); % Turns the image into a matrix

end

% The previous loop allows us to input the image either as a function argument

% (then CELLIMAGE is a string, for instance cell.jpg),

% or the image matrix (then CELLIMAGE is a matrix).

if nargin==1 % In case the dimensions of the matrix were not input as arguments.

[HEIGHT WIDTH]=size(MATRIX);

end

if nargin~=1 & nargin~=3

error (' Number of expected input arguments: 1 or 3');

end

WIND=MATRIX.*fspecial('gaussian',[HEIGHT WIDTH],0.2*min(HEIGHT,WIDTH)); % Windows the image with a gaussian filter

WIND=grayscale(WIND); % Rescales the entries between 0 and 255.

FOURTRANS=fft2(WIND); % Calculates the 2D Fourier transform

ABSOLU=abs(FOURTRANS); % Calculates the modulus of the Fourier coefficients

POWSPECDEN=ABSOLU.^2;

LOGA=log(POWSPECDEN);

FINAL=fftshift(LOGA); % Brings the main frequency peak in the center od the image

FINAL=grayscale(FINAL);

PSD display: disppsd.m(download)

function disppsd(MATRIX,FINAL)

% Displays the cellular pattern and the PSD.

% EXAMPLES:

% >> disppsd('cell.jpg')

% where cell.jpg is the cellular pattern image.

% OR

% >> disppsd(MATRIX)

% where MATRIX is the double-type matrix of

% the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg'));

% OR

% >> disppsd(MATRIX,FINAL)

% where MATRIX is the same as above and FINAL is the result of

% the PSD computation, i.e. FINAL=grayscale(calcpsd(MATRIX));

if isstr(MATRIX)

MATRIX=double(imread(MATRIX)); % Turns the image into a matrix

end

% The previous loop allows us to input the image either as a function argument

% (then CELLIMAGE is a string, for instance cell.jpg),

% or the image matrix (then CELLIMAGE is a matrix).

if nargin==1

FINAL=grayscale(calcpsd(MATRIX));

end

figure;colormap(gray);imshow(grayscale(MATRIX));image(grayscale(MATRIX));title('Disppsd -- Cellular Pattern');axis on;zoom on; % Imshow>>Image gives better results than Image alone. I don't know why...

figure;colormap(gray);imshow(FINAL);image(FINAL);title('Disppsd -- Power Spectral Density');axis on;zoom on;

Angular scanning: angplot.m(download)

function angplot(FINAL,HEIGHT,WIDTH,Xc,Yc)

% Puts intensity as a function of angle at defined radius.

% The result is a plot of the intensity along a circle

% centered on the crossing of the PSD-legs, at the radius

% the user is invited to input.

% The computation processes the average intensity on a window

% whose size can also be specified by the user.

% You may want to display the PSD before, to help you visualize

% and choose the appropriate radius.

% The main purpose of this program is to know the position of the PSD-legs.

% Once you know it, you can filter the PSD or run a radial scanning (radplot)

% along one of the legs.

% EXAMPLES:

% >> angplot(FINAL)

% where FINAL is the PSD matrix.

% You can also specify the HEIGHT, the WIDTH, and the center coordinates Xc and Yc

% of your image. If you do so, type:

% >> angplot(FINAL,HEIGHT,WIDTH,Xc,Yc)

% Enter either 1 argument (i.e. FINAL), or 5.

if nargin==1

[HEIGHT WIDTH]=size(FINAL);

[Xc Yc]=center(FINAL);

end

if HEIGHT~=WIDTH

disp(' WARNING: Your image isn''t square.')

disp(' The result this program is going to return are WRONG...')

end

if nargin~=1 & nargin~=5

error (' Number of expected input arguments: 1 or 5');

end

pi=3.14159265;

T=zeros(1,1024); % Intensity storage vector. Entries run through a circle at radius R0.

BUF=[]; % Temporary storage array

m=0;

while m~=3

if m==0 | m==1

iwid=input('\n Enter pixel window half-width for smoothing : ');

sprintf(' The averaging window size is now %dx%d',2*iwid+1,2*iwid+1)

sizeofwindow=(2*iwid+1)^2;

end

if m==0 | m==2

R0=input(' Enter the radius (in pixels) of the circle \n you want to run through : ');

end

h=waitbar(0,'Angular scanning in progress...not that easy!');

for I=1:1024

waitbar(I/1024,h)

TH=pi/512*(I-512); % Angle of study, from -pi to +pi

CT=cos(TH);

ST=sin(TH);

X=Xc-R0*ST; % Real coordinates of study

Y=Yc+R0*CT;

M=fix(X); % Top left closest pixel from the real point of study

N=fix(Y);

% CAUTION: the origin is at (1,1), in the top left corner of the image.

% That's why the expressions for X and Y are not those expected...

% Thus, the image is described in the trigonometric way,

% beginning on the right-hand side of the center.

if 1<=M & M<=HEIGHT & 1<=N & N<=WIDTH % Stay away from edges

for II=M:M+1 % Calculation around 4 pixels surrounding the real point

for J=N:N+1

LOCALSUM=0;

for K=-iwid:iwid % Average of the window intensity

for L=-iwid:iwid

if 1<=II+K & II+K<=HEIGHT & 1<=J+L & J+L<=WIDTH % Prevents the window from overlapping out of the image

LOCALSUM=LOCALSUM+FINAL(II+K,J+L)/sizeofwindow;

end

BUF(II,J)=LOCALSUM;

end

% Now we'got 4 values : one on each pixel surrounding

% the real point of study.

% To synthetize those 4 values into 1, we will average them by

% weighing them with the area of the opposite little square

% around the real point (the closer the pixel is from the point,

% the more important it is in the average).

P=X-M; % Lengths of the 4 squares

Q=Y-N;

PP=1-P;

QQ=1-Q;

T(I)=PP*QQ*BUF(M,N)+P*QQ*BUF(M+1,N)+Q*PP*BUF(M,N+1)+P*Q*BUF(M+1,N+1);

end

close(h)

figure;

I=1:1024;

ang=180/512*(I-512);

plot(ang,T)

axis on

xlabel('Angle in degrees, 0 is on the right hand-side of the center')

ylabel('Average intensity')

title(sprintf('Angplot -- Intensity along the circle of radius %d pixels -- Averaging window size : %dx%d',R0,2*iwid+1,2*iwid+1))

zoom on

m=menu('WHAT WOULD YOU LIKE TO DO NOW?','Change window averaging size','Change the scanning radius','Quit or return to main menu');

end

return

Radial scanning: radplot.m(download)

function ANGL=radplot(FINAL,HEIGHT,WIDTH,Xc,Yc)

% Puts intensity as a function of radius at defined angle.

% The result is a plot of the intensity along a leg of the PSD,

% centered on the crossing of the PSD-legs, at the angle

% the user is invited to input (this angle is set as an output,

% in case it is needed afterwards).

% The computation processes the average intensity on a window

% whose size can also be specified by the user.

% You may want to run the program angplot before, to know the accurate

% position of the PSD-legs.

% EXAMPLES:

% >> radplot(FINAL)

% where FINAL is the PSD matrix.

% You can also specify the HEIGHT, the WIDTH, and the center coordinates Xc and Yc

% of your image. If you do so, type:

% >> angplot(FINAL,HEIGHT,WIDTH,Xc,Yc)

% Enter either 1 argument (i.e. FINAL), or 5.

% To get the angle of the scanned leg and store it in THETA, type:

% >> THETA=radplot(...);

if nargin==1

[HEIGHT WIDTH]=size(FINAL);

[Xc Yc]=center(FINAL);

end

if nargin~=1 & nargin~=5

error (' Number of expected input arguments: 1 or 5');

end

pi=3.14159265;

T=zeros(1,2049); % Intensity storage vector. Entries run along a line at angle TH0.

BUF=[]; % Temporary storage array

m=0;

while m~=3

if m==0 | m==1

iwid=input('\n Enter pixel window half-width for smoothing : ');

sprintf(' The window averaging size is now %dx%d',2*iwid+1,2*iwid+1)

sizeofwindow=(2*iwid+1)^2;

end

if m==0 | m==2

TH0=input(' Enter the angle (in degrees) of the leg \n you want to run along : ');

ANGL=TH0; % Just for the convenience. ANGL for the output, TH0 in the program.

end

h=waitbar(0,'Radial scanning in progress...quite tough!');

for I=-1024:1024

waitbar((I+1024)/2048,h)

R=sqrt((HEIGHT/2)^2+(WIDTH/2)^2)/1024*I; % Radius of study

CT=cos(TH0*2*pi/360);

ST=sin(TH0*2*pi/360);

X=Xc-R*ST; % Real coordinates of study

Y=Yc+R*CT;

M=fix(X); % Top left closest pixel from the real point of study

N=fix(Y);

if 1<=M & M<=HEIGHT & 1<=N & N<=WIDTH % Stay away from edges

for II=M:M+1 % Calculation around 4 pixels surrounding the real point

for J=N:N+1

LOCALSUM=0;

for K=-iwid:iwid % Average of the window intensity

for L=-iwid:iwid

if 1<=II+K & II+K<=HEIGHT & 1<=J+L & J+L<=WIDTH % Prevents the window from overlapping out of the image

LOCALSUM=LOCALSUM+FINAL(II+K,J+L)/sizeofwindow;

end

BUF(II,J)=LOCALSUM;

end

% Now we'got 4 values : one on each pixel surrounding

% the real point of study.

% To synthetize those 4 values into 1, we will average them by

% weighing them with the area of the opposite little square

% around the real point (the closer the pixel is from the point,

% the more important it is in the average).

P=X-M; % Lengths of the 4 squares

Q=Y-N;

PP=1-P;

QQ=1-Q;

T(I+1024)=PP*QQ*BUF(M,N)+P*QQ*BUF(M+1,N)+Q*PP*BUF(M,N+1)+P*Q*BUF(M+1,N+1);

end

close(h)

figure;

I=-1024:1024; rad=(sqrt((HEIGHT/2)^2+(WIDTH/2)^2)/1024)*I;

plot(rad,T)

axis on

xlabel('Radius in pixels from the center')

ylabel('Average intensity')

title(sprintf('Radplot -- Intensity along the line of angle %d degrees -- Averaging window size : %dx%d',TH0,2*iwid+1,2*iwid+1))

zoom on

m=menu('WHAT WOULD YOU LIKE TO DO NOW?','Change window averaging size','Change the scanning angle','Quit or return to main menu');

end

return

Ultimate cell size computation: peak.m(download)

function RESULT=peak(DIMENSION,ANGL)

% Returns the sequence of the harmonics for the horizontal cell width.

% EXAMPLE:

% >> peak(DIMENSION,ANGL)

% where DIMENSION is the dimension of the square image

% and ANGL the angle of the PSD leg you've scanned.

m=2;

while m==2

BUF=input('\n Enter the sequence of the coordinates of the peaks, but not the DC! \n General syntax: [Abscissa#1 Ordinate#1;Abscissa#2 Ordinate#2; ... ] \n Example: [5 97;11 45;15 32] \n Now your turn: ');

X=BUF(:,1)';

Y=BUF(:,2)';

SIZEPIXELPERP=DIMENSION./X; % Cell size perpendiculary to the cells.

WEIGHT=(Y-max(Y))*(100-10)/(max(Y)-min(Y))+100; % Rescales the energy of the sequence: energy of lowest peak = 10, energy of highest peak = 100.

SIZEPIXELHORI=SIZEPIXELPERP/cos(ANGL*2*pi/360); % Projects the perpendicular cell width to the horizontal direction.

disp(' ');

disp(' Sequence of the harmonics for the horizontal cell size in pixels:');

disp(SIZEPIXELHORI);

m=menu('Does it make sense?','Yes, convert into centimeters','No, change values');

if m==1

disp(' ');

disp(' Relation between pixels and cm -- Use Photoshop ruler!');

c=input(' Enter a distance in cm: ');

p=input(' Enter the corresponding number of pixels: ');

SIZECMHORI=c*SIZEPIXELHORI/p;

RESULT=[SIZECMHORI' WEIGHT'];

disp(' ');

disp('**************************************************************');

disp(' ');

disp(' Here''s the result:')

disp(' 1st column: horizontal cell size in cm, and its harmonics.');

disp(' 2nd column: relative energy of each harmonic (linear fit):');

disp(' Energy of lowest peak = 10, energy of highest peak = 100');

disp(' ');

disp(RESULT);

disp(' ');

disp('**************************************************************');

disp(' ');

end

return

Detection of the center of the PSD: center.m(download)

function [Xc,Yc]=center(M)

% >> [Xc,Yc]=center(MATRIX) returns as Xc and Yc the row

% and the column of the pixel of highest

% intensity in the matrix MATRIX.

% Usefull to find the center of a PSD image!

[m,I]=max(M);

[mm,Yc]=max(m);

Xc=I(Yc);

Picture grayscaling: grayscale.m(download)

function GRAYSCALED=grayscale(MATRIX)

% Returns a matrix whose entries have been lineary

% stretched between the initial minimum and maximum

% entries, in order to occupy the whole gray scale,

% from 0 to 255.

% EXAMPLE:

% >> GRAYSCALED=grayscale(MATRIX)

M=max(max(MATRIX));

m=min(min(MATRIX));

GRAYSCALED=255/(M-m)*(MATRIX-m);

cellsizer.m.........	main program
filtr.m................	filtering program
matlabsobel.m...	Matlab-based edge enhancement program
persobel.m........	home-made edge enhancement program
calcpsd.m.........	PSD computation
disppsd.m.........	PSD display
angplot.m..........	angular scanning
radplot.m..........	radial scanning
peak.m..............	ultimate cell size computation
center.m............	detection of the center of the PSD
grayscale.m.......	picture grayscaling