Matlab: map an image with lookup table

Suppose we have a lookup table (lut.txt) look like (two columns):

0     0
1     0.0005
2     0.0006
3     0.002
.
.
254   0.981
255   1

Only save the 2nd column to LUT0.txt:

0
0.0005
.
.
0.981
1

matlab code:

load LUT0.txt;
img_gamma = LUT0(im2uint8(InputImage)+1);

solve equation in MATLAB

>> syms x y d1 d2 R1 R2 mag1 mag2 mag3;

>> S = solve('x^2+y^2=R1^2', '(d2-x)^2+y^2=R1^2*mag1/mag2','x^2+(d1-y)^2=R1^2*mag1/mag3','x','y','R1');

>> S.x
>> S.y

Request user input in MATLAB

Wait for typing from keyboard

var = input('what is the value:\n');
filename = input ('what is the name?', 's');
 
if isempty(filename)
    filename = 'default.dat'
end

Matlab: Sound generation

Pure tone: unmodulated sinusoid, matlab function:

% Ft: frequency of the tone
% Fs: sample rate
% Td: last seconds
% Output is a row vector s containing the sampled point
% t is an optional and is the time axis associated with s

function [s, t] = tonegen(Ft, Fs, Td)
t = [1:Fs*Td]/Fs;
s = sin (2 * pi * t * Ft);

% example:
y = tongen(440, 16000, 5);
soundsc(y, 16000);

Generate sound with different frequency:

% frc: time-varying frequency
% Fs: sample rate

function  [snd]=freqgen(frc,  Fs)
th=0;
fr=frc*2*pi/Fs;
for  si=1:length(fr)
 th=th+fr(si);
 snd(si)=sin(th);
 th=unwrap(th);
end

% example
>> freq=[440*(1+zeros(1,1000)),415.2*(1+zeros(1,1000)),392*(1+zeros(1,1000))];
>> music=freqgen(freq,  8000);
>> soundsc(music,  8000);

Audio in Matlab: record, play, and display

Recording using ‘wavrecord()’:

>> speech = wavrecord(16000, 8000, 1, 'double');
>> sound(speech) % play sound with default sampling rate 8192Hz)

This records 16000 samples with a sample rate of 8kHz in mono (‘1’). (Windows only).

recording using ‘audiorecorder()’:

>> x = audiorecorder(16000, 16, 1);
>> record(x);
% recording ...
>> stop(x); % other commands: pause(), resume()
>> play(x);
>> speech2 = getaudiodata(x, 'double'); %convert to normal vector

x: a handle of audio recorder object.
16000: Hz, sampling rate (common rates are:8k,11025,22050,44000)
16: sample precision in bits
1: number of channels

The max range values in 16-bit format is -32768~32767. It is scaled to –1.0~1.0 for double precision. If for some reason it is out of range, use following to play instead of ‘sound(speech)’:

>> soundsc(speech, 8000); % or
>> sound(speech / max(abs(speech)), 8000);

Display in time-domain:

>> plot(speech); % or
>> plot([1:size(speech)]/8000, speech); % time in second in x-axis

Frequency response of Laplacian filter

zero freq response shifted to (32,32), Matlab code:

clear
X=[-1 0 -1;-1 4 -1;-1 0 -1];
Y=fft2(X,64,64);
Y=fftshift(Y);
Y=abs(Y);
figure(1)
[x,y]=meshgrid(1:1:64,1:1:64);
mesh(x,y,Y)
title('frequency response of 3x3 Laplacian filter')
xlabel('k_1 +32')
ylabel('k_2 +32')
figure(2)
contour(Y)
title('frequency response of 3x3 Laplacian filter')
xlabel('k_1 +32')
ylabel('k_2 +32')

MATLAB: bwselect, url

blob = bwselect(bw);
imshow(blob);

the function displays bw and waits for you to click on the image with the mouse. You can click on one or multiple points. To finish your selection, double-click on the last point, or just press RETURN. bwselect then returns a binary image containing only the objects ("blobs") that you clicked on. You can also chose one position (10, 23) with:

blob = bwselect(bw, 10, 23);
clf % clear figure;
imshow(blob)

Read image from Internet link:

url = 'http://abc.com/images/us.png';
bw = imread(url);
imshow(bw)

MATLAB: deal with .avi file

Read the video file and get all the information:

mov = aviread('carslave.avi');
movinfo = aviinfo('carslave.avi')
noframe = movinfo.NumFrames;
movie(mov)

Get the first frame and change to image file:

frame1 = aviread('carslave', 1);
img = frame2im(frame1);
imshow(img);

frame2im only deals with uint8 format. It doesn't work with some avi files with uint16 data. Mark here and try to find solution later.

Get the image files and play in movie:

for i = 1:100
...
[X, map] = gray2ind(img, 256); %rgb2ind
mov(i) = im2frame(X, map);
end
movie(mov);

Create avi file:

%aviobj = avifile('test.avi');
aviobj = avifile('test.avi','fps',25,'COMPRESSION','None');
for i = 1:100
...
[X, map] = gray2ind(img, 256);
rgbimg = ind2rgb(X, map);
aviobj = addframe(aviobj, rgbimg);
end
aviobj = close(aviobj);

Check this videoIO toolbox

MATLAB: is*

The following 'is-' are useful in image processing,

if isrgb(img)
isgray, isbw, isrgb, isind, isa

check with is* in matlab

Basic image compression technologies

1 . Block Truncation Coding (link)

1. divide image into 4x4 or 8x8 blocks
2. Calculate the average gray level of the block
3. use a threshold to get binary block b and mu, ml (greater and less then th)
4. The encoder writes mu, ml and b to a file
5. In the decoder, an image block is reconstructed by replacing the 1's with mu and the 0's by ml.

2. Gaussian Pyramids (link): Each level in the pyramid is 1/4 of the size of previous level.

1. low-pass filter the image,
2. down-sampling image,
3. up-sampling image,
4. low-pass filter again.

3. Singular Value Decomposition (link). A = USV^T, U and V are orthogonal matrices (left/right singular vectors), and S is a diagonal matrix (contains the singular values). The size of U S V is m \times m, mxn, nxn, respectively.
This method approximated the mxn image by using far fewer entries than in original matrix. By using the rank of a matrix, we remove the redundant information. Matlab code:

dvalue= im2double(imread('map.bmp'));
[u,s,v] = svds(dvalue, singvals); %ex. singvals=20
im = (u * s * transpose(v));

where 'svds(A,k)' computes the k largest singular values and associated singular vectors of matrix A.

4. DCT: former post on DCT and JPEG compression

MATLAB: numel, linspace, and nan

number of element

temp = [1 2; 3 4];
numel(temp)

ans =

4

linspace: generate linear spaced vectors:

linspace(1, 6, 9)

ans =

1.0000 1.6250 2.2500 2.8750 3.5000 4.1250 4.7500 5.3750 6.0000

create a matrix of non-number (not a number) and we can fill in it with number later:

nan(4, 3)

ans =

NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN

JPEG compression

The basic jpeg compression steps (only think of gray level image):

1. divide the image to 8x8 matrices
2. use DCT to each 8x8 matrix
3. sum up all the 8x8 matrices and get a matrix sum_mat
   
4. sort the elements of 8x8 matrix(sum_mat) from the highest to the smallest, and get the indices list
5. choose (some) highest coefficients: the ratio: (some high energy)/(total energy) is good enough (or pick the first N highest coefficients)
   
6. save the index list of the high energy (coefficient). In each 8x8 matrix, only save these coefficients.
7. build 8x8 matrices with zero and store the right coefficient in the right place with the index list.
   
8. use IDCT to each matrix
9. 8x8 matrices to big image size matrix. (restored_img)
   
10. Calculate SNR: SNR = (sum(original_img(:).^2)) / (sum(original_img(:)-restored_img(:)).^2)

Also look at the book digital image processing 2nd Ed. page 498. Some relevant matlab code can be found at former post Use DCT for image compression

Matlab code:Use DCT for image compression

Just show the idea for gray level image. RGB image is in similar way. It takes the following steps:

1. DCT conversion and create a 1D matrix
2. Sort elements of the matrix from the highest to the smallest and get the indices list
3. Choose the highest N elements
4. Inverse DCT

Matlab code:

dvalue=imread('map.bmp');
dvalue = double(dvalue)/255;
dvalue = rgb2gray(dvalue);
img_dct=dct2(dvalue);
img_pow=(img_dct).^2;
img_pow=img_pow(:);
[B,index]=sort(img_pow);%no zig-zag
B=flipud(B);
index=flipud(index);
compressed_dct=zeros(size(dvalue));
coeff = 20000;% maybe change the value
for k=1:coeff
compressed_dct(index(k))=img_dct(index(k));
end
im=idct2(compressed_dct);
imwrite(im, 'maptemp2.bmp')

MATLAB: figure tips

Always forgot which is x or y, confusing with 'img(y, x)--(row, column)'

figure(2)
[x, y, button] = ginput(2);
plot(x, y, 'r*');

make a figure current and make it ready for operation

figure(5) %'h=figure(5)': h=5
set(0, 'CurrentFigure', 5); % make figure current

Change the name (title) of the figure

h = figure('Name','my window','NumberTitle','off'); %h=1,2 ..N
set(0, 'CurrentFigure', h);

MATLAB: save poistions with mouse-click (ginput)

pt_cam = [];
im1 = imread( 'B7_3.bmp' );
figure(1);imshow(im1);hold on;
im2 = imread( 'map_B7_1.bmp' );
figure(2);imshow(im2);hold on;

figure(1);
[xc,yc,Button] = ginput(1);
plot(xc,yc,'r.','Markersize',15)
figure(2);
[xm,ym,Button] = ginput(1);
plot(xm,ym,'.','Markersize',15)

while Button~=3
pt = [xc yc];
pt_cam = [pt_cam; pt];

figure(1);
[xc,yc,Button] = ginput(1);
plot(xc,yc,'r.','Markersize',15)

if Button==3;
break;
end

figure(2);
[xm,ym,Button] = ginput(1);
plot(xm,ym,'.','Markersize',15)
end

MATLAB/Simulink works with Windows and TI TMS320DM642 on video processing

Windows platform:

1. Input video from an AVI file. We recommend 'image signal' to be 'One multidimensional signal'.
2. Use the first frame as the background image
3. embed MATLAB function: drag Simulink-User/Defined Functions/Embeded MATLAB Function, double-click to edit.
   

4. function y = fcn(u)
   y = false;
   k = 0;
   for i = 158:161
   for j = 220:223
   if u(i,j) > 0
   k = k+ 1;
   end
   end
   end
   if k>15
   y = true;
   end

5. embed C function: We can either write an S-function using the templates and examples or use the legacy code tool (lct). The legacy code tool is probably the best way to incorporate c-code into our model, the following m-file script shows how to incorporate c-code with the legacy code tool. (lct.m)

6. %% Initialize the structure necessary for using the legacy code tool
   def = legacy_code('initialize');
   
   %% Fill in the appropriate structure fields
   % We define the source files, header files, libraries, and associated paths
   % first. I am assuming the files are in the same directory as where this
   % script will run.
   def.SourceFiles = {'my_Fun.c'};
   def.HeaderFiles = {'my_fun.h'};
   
   % define the name for your s-function
   def.SFunctionName = 'my_Fun_sfun';
   
   % define the simulink equivalent function spec to the c-code. In our model,
   % the input and output were both boolean, so I changed the function spec to
   % reflect that. Also, by using u1[][], I am assuming a matrix of unknown
   % dimensions is passed to the function. This is the simplest way to set it
   % up.
   def.OutputFcnSpec = 'boolean y1 = my_Fun( boolean u1[][] )';
   
   %% Auto generate an S-function / compile it / create a block
   legacy_code('sfcn_cmex_generate', def);
   legacy_code('compile', def);
   legacy_code('slblock_generate', def);

7. Use 'Video Viewer' to check the output video.
8. To add some operations before Simulink codes running, right-click on blank area and choose 'Model Properties'/Callbacks. Now you can add something on 'PreLoadFun*'. We use the following m-file script to load an image and choose some positions with mouse-click:

9. f = figure;
   imshow('1st.bmp');
   title('Right Click to stop input');
   hold on;
   
   % Input features with mouse
   path = [];
   i = 0;
   [x,y,Button] = ginput(1);
   while Button~=3
   pt = [x y];
   path = [path; pt];
   plot(x,y,'g.','Markersize',15)
   if (i~=0)
   plot([x,temp(1)], [y, temp(2)], 'g-');
   end
   i = 1;
   % Give world coordinates of selected feature
   [x,y,Button] = ginput(1);
   temp = pt;
   end
   close(f);

10. All the mouse-click positions are transferred to 'Draw Polygons' to show the region.

TI TMS320DM642 platform:

1. MATLAB/Simulink version:2007b. DM642 version: version 3 with 720MHz, TVP5146 decoder. Code Composer Studio (CCS) version: 3.3. 2007b only works with CCS3.3. 2007a and older versions only work with CCS3.0,3.1 and they don't support CCS3.3.
2. Choose 'DM642EVM V3' from 'Target for TI C6000/C6000 Target Preferences', double check the CPU clock (720MHz)
3. Video Capture blockset: Decoder type:TVP5146, Input port:Port 0, NTSC. 'Analog video input', We am using 'S-video' and the Simulink demo uses Composite. 'Data order': 'Column major'.The text characters are flipped if choosing Row. The input source is Camera/DVD player
   
4. Video Display: same 'data order'.
5. MATLAB/Simulink Demo: 'Video Stabilization' and 'Edge Detection' are good demos to check your CCS , DM642 board and video in/out system. Both demos need 'fixed-point' blockset. 'Surveillance' demo does NOT use external camera input or video output. So nothing will show up on our monitor.
   
6. Without running anything on DSP, the output display should be color bars.
   
7. You can use "Concatenate" block to implement a quad-display for output.
8. the data type of the substraction of current image from background image is an int16. This was set to uint8 in Windows version and was causing underflow
9. Open CCS3.3 and make sure the board is connected.
   
10. Use 'incremental build' to build the project from CCS, load to DM642 and run.
11. The following message appear on MATLAB command window during the project building:

12. >>
    Target Preferences block specifies that the board named 'DM642 EVM XDS560 Emulator' will be used to run generated code.
    However, since only board named 'Blackhawk USB560m - DM642 EVM Emulator' is found in your system, that board will be used.
    ### Connecting to Code Composer Studio(tm) ...
    ### Generating code into build directory: C:\Malab_Eva\lq_TI_ccslink
    Embedded MATLAB parsing for model "lq_TI"...Done
    Embedded MATLAB code generation for model "lq_TI"....Done
    Embedded MATLAB compilation for model "lq_TI"... 1 file(s) copied.
    Done
    Embedded MATLAB parsing for model "lq_TI"...Done
    Embedded MATLAB code generation for model "lq_TI"....Done
    ### Invoking Target Language Compiler on lq_TI.rtw
    tlc
    -r
    C:\Malab_Eva\lq_TI_ccslink\lq_TI.rtw
    C:\Program Files\MATLAB\R2007b\toolbox\rtw\targets\ccslink\ccslink\ccslink_grt.tlc
    -OC:\Malab_Eva\lq_TI_ccslink
    -IC:\Program Files\MATLAB\R2007b\toolbox\rtw\targets\ccslink\ccslink
    -IC:\Program Files\MATLAB\R2007b\toolbox\vipblks\vipmex\tlc_c
    -IC:\Program Files\MATLAB\R2007b\toolbox\rtw\targets\tic6000\tic6000blks\tlc_c
    -IC:\Program Files\MATLAB\R2007b\toolbox\dspblks\dspmex\tlc_c
    -IC:\Program Files\MATLAB\R2007b\toolbox\shared\etargets\etargets\tlc_c
    -IC:\Malab_Eva\lq_TI_ccslink\tlc
    -IC:\Program Files\MATLAB\R2007b\rtw\c\tlc\mw
    -IC:\Program Files\MATLAB\R2007b\rtw\c\tlc\lib
    -IC:\Program Files\MATLAB\R2007b\rtw\c\tlc\blocks
    -IC:\Program Files\MATLAB\R2007b\rtw\c\tlc\fixpt
    -IC:\Program Files\MATLAB\R2007b\stateflow\c\tlc
    -aEnforceIntegerDowncast=1
    -aFoldNonRolledExpr=1
    -aInlineInvariantSignals=1
    -aInlineParameters=0
    -aLocalBlockOutputs=1
    -aRollThreshold=5
    -aZeroInternalMemoryAtStartup=1
    -aZeroExternalMemoryAtStartup=1
    -aInitFltsAndDblsToZero=1
    -aGenerateReport=0
    -aGenCodeOnly=1
    -aRTWVerbose=1
    -aIncludeHyperlinkInReport=0
    -aLaunchReport=0
    -aGenerateTraceInfo=0
    -aForceParamTrailComments=0
    -aGenerateComments=1
    -aIgnoreCustomStorageClasses=1
    -aIncHierarchyInIds=0
    -aMaxRTWIdLen=31
    -aShowEliminatedStatements=0
    -aIncDataTypeInIds=0
    -aInsertBlockDesc=0
    -aSimulinkBlockComments=1
    -aInlinedPrmAccess="Literals"
    -aTargetFcnLib="ansi_tfl_tmw_new.mat"
    -aIsPILTarget=0
    -aLogVarNameModifier="rt_"
    -aGenerateFullHeader=1
    -aexportIDEObj=1
    -aideObjName="CCS_Obj"
    -aProfileGenCode=1
    -aInlineDSPBlks=1
    -abuildAction="Build_and_execute"
    -acompilerOptionsStr="-o2 --mem_model:data=far"
    -alinkerOptionsStr=""
    -asystemStackSize=8192
    -aoverrunNotificationMethod="None"
    -aoverrunNotificationFcn="myfunction"
    -aGenerateTraceInfo=0
    -aProfileNumSamples=100
    -asystemStackSize=8192
    -p10000
    
    
    ### Loading TLC function libraries
    
    ......
    ### Initial pass through model to cache user defined code
    ....
    ### Caching model source code
    .........................................................
    ### Writing header file lq_TI_types.h
    ### Writing header file lq_TI.h
    ### Writing header file lq_TI_private.h
    ### Writing source file lq_TI.c
    .
    ### Writing header file rtmodel.h
    ### Writing source file lq_TI_data.c
    ### Writing header file rt_nonfinite.h
    ### Writing source file rt_nonfinite.c.
    
    ### Writing source file lq_TI_main.c
    ### Writing header file lq_TI_main.h
    ### Writing header file MW_c6xxx_csl.h
    .
    ### Writing source file MW_c6xxx_csl.c
    ### TLC code generation complete.
    ### Creating project ...
    ### Building Code Composer Studio(tm) project...
    ### Build complete
    ### Downloading COFF file
    ### Downloaded: lq_TI.out

13. CCS starts running the code on DM642

Problems:

• In embedded C-function, how to make the input parameter to be an image array
• Sometimes the program is extremely slow on Windows platform (7-8 frames/second).
  
• How to check frame rate on DM642 platform.
• How to make the debugging process more convenient, such as printing the variable values and image pixel values. Update from ICIP2007: highlight the data line and we can add 'mplayer' on it to show the pixel values. For variable values, we may add a blockset (scope?) to check the value in running time.
  
• When will MATLAB/Simulink support DaVinci.
  

Debug command in Matlab

This can be used after 'M-x matlab-shell' in emacs to debug the matlab code

dbstop laser.m 12
run laser
dbstep (dbstep 13)
dbstep in (go inside function, dbstep out)
dbclear % clear breakpoints
dbcont % resume execution
dbquit
dbstatus % list all breakpoints

Save date to text file with 'fprintf'

Save the data to several columns (x, y, z)

x = 1:10;
y = 11:20;
fid = fopen('temp2.txt','wt');
fprintf(fid, '%d %d \n', x, y);
fclose(fid);

The result look like:

1 11
2 12
3 13
4 14
5 15
6 16
7 17
8 18
9 19
10 20

Similiar method:

x = 1:10;
y = 11:20;
s(1,:) = x;
s(2,:) = y;
fid = fopen('temp2.txt','wt');
fprintf(fid, '%d %d \n', s);
fclose(fid);

inverse z transform

syms z n;
fz = 5+2*z^(-1)+4*z^(-3);
iztrans(fz)


/>ans =
/>5*charfcn[0](n)+2*charfcn[1](n)+4*charfcn[3](n)

where charfcn[1](n) is the delta function: [tex]\delta[n-1] [/tex]

imdistline, imline, impoint, and imrect

Check the following several functions. I used my own scripts before.

imdistline, imline, impoint,
and imrect