1 | %% Copyright (c) 2001, NEC Research Institute Inc. All rights reserved.
|
---|
2 | %%
|
---|
3 | %% Permission to use, copy, modify, and distribute this software and its
|
---|
4 | %% associated documentation for non-commercial purposes is hereby
|
---|
5 | %% granted, provided that the above copyright notice appears in all
|
---|
6 | %% copies, derivative works or modified versions of the software and any
|
---|
7 | %% portions thereof, and that both the copyright notice and this
|
---|
8 | %% permission notice appear in the documentation. NEC Research Institute
|
---|
9 | %% Inc. shall be given a copy of any such derivative work or modified
|
---|
10 | %% version of the software and NEC Research Institute Inc. and its
|
---|
11 | %% affiliated companies (collectively referred to as NECI) shall be
|
---|
12 | %% granted permission to use, copy, modify and distribute the software
|
---|
13 | %% for internal use and research. The name of NEC Research Institute
|
---|
14 | %% Inc. and its affiliated companies shall not be used in advertising or
|
---|
15 | %% publicity related to the distribution of the software, without the
|
---|
16 | %% prior written consent of NECI. All copies, derivative works or
|
---|
17 | %% modified versions of the software shall be exported or reexported in
|
---|
18 | %% accordance with applicable laws and regulations relating to export
|
---|
19 | %% control. This software is experimental. NECI does not make any
|
---|
20 | %% representations regarding the suitability of this software for any
|
---|
21 | %% purpose and NECI will not support the software. THE SOFTWARE IS
|
---|
22 | %% PROVIDED AS IS. NECI DOES NOT MAKE ANY WARRANTIES EITHER EXPRESS OR
|
---|
23 | %% IMPLIED WITH REGARD TO THE SOFTWARE. NECI ALSO DISCLAIMS ANY WARRANTY
|
---|
24 | %% THAT THE SOFTWARE IS FREE OF INFRINGEMENT OF ANY INTELLECTUAL PROPERTY
|
---|
25 | %% RIGHTS OF OTHERS. NO OTHER LICENSE EXPRESS OR IMPLIED IS HEREBY
|
---|
26 | %% GRANTED. NECI SHALL NOT BE LIABLE FOR ANY DAMAGES, INCLUDING GENERAL,
|
---|
27 | %% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, ARISING OUT OF THE USE
|
---|
28 | %% OR INABILITY TO USE THE SOFTWARE.
|
---|
29 |
|
---|
30 |
|
---|
31 | % This README file will explain a little bit about my code for linear fitting
|
---|
32 | % with missing data. This code addresses the problem of fitting a low rank
|
---|
33 | % approximation to a matrix in which some of the elements are unknown. That is,
|
---|
34 | % this algorithm attacks the problem that is easily solved by SVD when there
|
---|
35 | % is no missing data. This directory also contains code that is specialized to
|
---|
36 | % solving structure from motion (SFM) problems. A full description of the algorithms
|
---|
37 | % can be found in ``Linear Fitting with Missing Data for Structure-From-Motion,''
|
---|
38 | % Computer Vision and Image Understanding, 82:57-81, (2001), D. Jacobs.
|
---|
39 | % An earlier version of this work, with earlier code, was released in '97 and
|
---|
40 | % presented in CVPR '97. That code did not have the specializations for SFM,
|
---|
41 | % and did not work on motion sequences with translation (which requires a bit
|
---|
42 | % more than a generic low rank approximation).
|
---|
43 |
|
---|
44 | % This directory contains not just the code for my algorithm, but also the
|
---|
45 | % code used to generate all the tests described in the CVIU paper, and
|
---|
46 | % some additional tests as well. In particular, there is also code that
|
---|
47 | % implements an iterative method due to Wiberg, which is described in a paper
|
---|
48 | % by Shum, Ikeuchi and Reddy in PAMI, Sept. 1995. And there is code to
|
---|
49 | % generate synthetic motion sequences, occlusions and intensity images. I am
|
---|
50 | % only trying to introduce users to the central code in this file; if you want
|
---|
51 | % to try to use the code developed for my experiments you'll have to nose around
|
---|
52 | % these files.
|
---|
53 |
|
---|
54 | % Note that my earlier code (rankr.m) does have a lot of little hacks that are not
|
---|
55 | % well documented. For the most part these concern choosing column triples
|
---|
56 | % to estimate the nullspace, and methods for extending an initial stable solution
|
---|
57 | % when the stable solution only covers part of the matrix. This typically occurs
|
---|
58 | % when the amount of missing data is large.
|
---|
59 |
|
---|
60 | % COMPATIBILITY NOTE: This code has been upgraded to run in MATLAB V. As
|
---|
61 | % a consequence, it won't run in MATLAB IV. To get it to run in MATLAB IV
|
---|
62 | % you must:
|
---|
63 |
|
---|
64 | % -- Remove the call to "logical" in whole_invert.
|
---|
65 |
|
---|
66 | % On the other hand, I have not completely tested the port to MATLAB V.
|
---|
67 | % However, the main functions, shown below, have been tested as shown.
|
---|
68 |
|
---|
69 | % This file is also a Matlab M-file, and can be run to see a quick demo of
|
---|
70 | % some of the main functions.
|
---|
71 |
|
---|
72 | % Although I don't promise to support this code, or upgrade, any questions
|
---|
73 | % or comments can reach me at: djacobs@cs.umd.edu.
|
---|
74 |
|
---|
75 | % First, I'll demonstrate the generic code on a special case for which it is
|
---|
76 | % particularly good.
|
---|
77 |
|
---|
78 | M = 100.*(rand(6,3)*(rand(3,9)))
|
---|
79 | INC = [1,1,1,1,1,1,0,0,0; ...
|
---|
80 | 1,1,1,1,1,1,0,0,0; ...
|
---|
81 | 1,1,1,0,0,0,1,1,1; ...
|
---|
82 | 1,1,1,0,0,0,1,1,1; ...
|
---|
83 | 0,0,0,1,1,1,1,1,1; ...
|
---|
84 | 0,0,0,1,1,1,1,1,1]
|
---|
85 |
|
---|
86 | % We've just generated a small, 6x9 matrix that is rank 3.
|
---|
87 | % INC provides a pattern of missing data for this matrix, with
|
---|
88 | % 1 indicating the data is present, and 0 that it is missing.
|
---|
89 | % None of the algorithms will look at the data associated with
|
---|
90 | % 0s, except the code to evaluate the results. This matrix has
|
---|
91 | % the interesting property that there is no 4x4 subblock that is missing
|
---|
92 | % only one element, so the Tomasi and Kanade approach of "hallucinating"
|
---|
93 | % the values of missing elements will not work. To apply my algorithm
|
---|
94 | % to this data, one can type:
|
---|
95 |
|
---|
96 | disp('Our solution is:');
|
---|
97 | [e,a] = rankr(M,INC,3,1000)
|
---|
98 |
|
---|
99 | % This finds a rank 3 approximation to M, with the values in INC missing.
|
---|
100 | % The fourth parameter, 1000, says that 1000 triples of columns are the maximum
|
---|
101 | % number that will be used in estimating the nullspace; if we don't have enough
|
---|
102 | % information after that, we give up (see the paper if this doesn't make
|
---|
103 | % any sense). For this particular example, it would make more sense to just try
|
---|
104 | % all triples of columns, instead of a random sample, but for big matrices that
|
---|
105 | % won't work. However, because of this randomness, there's a small chance that the
|
---|
106 | % algorithm won't find the right answer.
|
---|
107 | % The value e that is returned is the sum of square differences between the
|
---|
108 | % elements that are present in the matrix, and the corresponding values in the
|
---|
109 | % approximating matrix. a is the approximating matrix.
|
---|
110 |
|
---|
111 | % This matrix also has too many missing elements for the iterative method of
|
---|
112 | % Shum et al to work. Here's how we call it:
|
---|
113 |
|
---|
114 | disp('The iterative method gives:');
|
---|
115 | [e,a] = shum(M,INC,3,0,100)
|
---|
116 |
|
---|
117 | % Here, e = -2 signals that the algorithm couldn't find a solution.
|
---|
118 | % The parameter 3 means find a rank 3 approximation, the parameter 0
|
---|
119 | % means pick a random matrix as the starting point (we could use a starting guess
|
---|
120 | % here), and 100 means that if the algorithm doesn't converge after 100 iterations,
|
---|
121 | % stop anyway.
|
---|
122 |
|
---|
123 | S = random_motion(4);
|
---|
124 | P = [rand(3,12); ones(1,12)];
|
---|
125 | % S contains random motions for 4 frames. P contains 12 random points.
|
---|
126 | % The ones are there so S*P gives the image points.
|
---|
127 | disp('M are image points from a random motion');
|
---|
128 | M = S*P
|
---|
129 | INC = [1,1,1,1,1,1,1,1,1,0,0,0;...
|
---|
130 | 1,1,1,1,1,1,1,1,1,0,0,0;...
|
---|
131 | 1,1,1,1,1,1,0,0,0,1,1,1;...
|
---|
132 | 1,1,1,1,1,1,0,0,0,1,1,1;...
|
---|
133 | 1,1,1,0,0,0,1,1,1,1,1,1;...
|
---|
134 | 1,1,1,0,0,0,1,1,1,1,1,1;...
|
---|
135 | 0,0,0,1,1,1,1,1,1,1,1,1;...
|
---|
136 | 0,0,0,1,1,1,1,1,1,1,1,1]
|
---|
137 | % This also represents a pattern of missing data that our algorithm can handle that
|
---|
138 | % is tricky for other approaches. See the paper for details.
|
---|
139 |
|
---|
140 | [e,a,s] = rankrsfm(M,INC)
|
---|
141 |
|
---|
142 |
|
---|
143 | % Now, let's try an example that all methods work on.
|
---|
144 | M = rand(20,3)*rand(3,20);
|
---|
145 | disp('We set M to be a rank 3, 20x20 matrix.');
|
---|
146 | INC = rand(20,20)>.25;
|
---|
147 | disp('We set INC to be a random matrix, with each element 1 with prob. .75.');
|
---|
148 |
|
---|
149 | disp('Our solution has error of:');
|
---|
150 | [e,a] = rankr(M,INC,3,1000);
|
---|
151 | disp(e);
|
---|
152 |
|
---|
153 | disp('If we use it as the starting point for an iterative method, we get error of:');
|
---|
154 | [e1,a1] = shum(M,INC,3,a,100);
|
---|
155 | disp(e1);
|
---|
156 |
|
---|
157 | disp('If we run the iterative method with a random starting point, we get error of:');
|
---|
158 | [e2,a2] = shum(M,INC,3,0,100);
|
---|
159 | disp(e2);
|
---|
160 |
|
---|
161 | % Although results are randomized, probably all three approaches gave a good result
|
---|
162 | % here.
|
---|
163 |
|
---|
164 | % There is also a lot of code to run experiments, generate synthetic test data, and
|
---|
165 | % compare methods. If you look at compare_motion, you'll get a starting point to
|
---|
166 | % the code I used to test the generic method for the CVPR '97 paper. compare_transmot
|
---|
167 | % contains code used to run experiments on the new SFM code.
|
---|
168 |
|
---|
169 | % Good luck
|
---|
170 | % -- David Jacobs
|
---|