source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/image3dstiching/Occlu/FindOccluPoint.m @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d
File size: 11.0 KB
Line 
1% *  This code was used in the following articles:
2% *  [1] Learning 3-D Scene Structure from a Single Still Image,
3% *      Ashutosh Saxena, Min Sun, Andrew Y. Ng,
4% *      In ICCV workshop on 3D Representation for Recognition (3dRR-07), 2007.
5% *      (best paper)
6% *  [2] 3-D Reconstruction from Sparse Views using Monocular Vision,
7% *      Ashutosh Saxena, Min Sun, Andrew Y. Ng,
8% *      In ICCV workshop on Virtual Representations and Modeling
9% *      of Large-scale environments (VRML), 2007.
10% *  [3] 3-D Depth Reconstruction from a Single Still Image,
11% *      Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng.
12% *      International Journal of Computer Vision (IJCV), Aug 2007.
13% *  [6] Learning Depth from Single Monocular Images,
14% *      Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng.
15% *      In Neural Information Processing Systems (NIPS) 18, 2005.
16% *
17% *  These articles are available at:
18% *  http://make3d.stanford.edu/publications
19% *
20% *  We request that you cite the papers [1], [3] and [6] in any of
21% *  your reports that uses this code.
22% *  Further, if you use the code in image3dstiching/ (multiple image version),
23% *  then please cite [2].
24% * 
25% *  If you use the code in third_party/, then PLEASE CITE and follow the
26% *  LICENSE OF THE CORRESPONDING THIRD PARTY CODE.
27% *
28% *  Finally, this code is for non-commercial use only.  For further
29% *  information and to obtain a copy of the license, see
30% *
31% *  http://make3d.stanford.edu/publications/code
32% *
33% *  Also, the software distributed under the License is distributed on an
34% * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
35% *  express or implied.   See the License for the specific language governing
36% *  permissions and limitations under the License.
37% *
38% */
39function [POriReprojM TargetID FieldID FieldOccluPix OcluedDist FaceSetIDRemained] = FindOccluPoint(defaultPara, ScaleImg, ITarget, IField, F, ...
40                TargetPointPix, FieldPointPix, INDTarget, INDField, ModelInfoTarget, ModelInfoField, Pair)
41
42% Function find the occlu point
43% And output a region on epipolar line
44
45% initial parameters
46FlagDisp = 1;
47EpipolarThre = 0.1;
48MaxFieldPointDist2Ray = 1;
49[dump DepthScale(1) DepthScale(2)] = size(ModelInfoField.Model.PlaneParaInfo.Position3DFited);
50
51NumTargetPointPix = size(TargetPointPix,2);
52AllFieldPoint3D = ModelInfoField.Model.PlaneParaInfo.Position3DFited;
53AllFieldPoint3D = AllFieldPoint3D(:,:);
54
55% change FieldPoint3D into Target Coordinate
56R_f_t = Pair.R';
57T_f_t = -Pair.R'*Pair.T;
58AllFieldPoint3D = R_f_t*AllFieldPoint3D + repmat( T_f_t, 1, size(AllFieldPoint3D,2)); % now in Target coordinate
59
60% Initialize Variables
61POriReprojM = [];
62TargetID = [];
63FieldID = [];
64FaceSetIDRemained = [];
65FieldOccluPix = [];
66OcluedDist = [];
67
68% Test Occlusion for each TargetPointPix point
69fprintf(' Finding Occlusion Points ......');
70FindOccluPointTime = tic;
71for IndexTargetPointPix = 1:NumTargetPointPix
72   
73%     if IndexTargetPointPix/100 == round(IndexTargetPointPix/100)
74%            IndexTargetPointPix
75%     end
76        % Initalize temp variables
77        INDFieldCandidate = [];
78
79    % 1) Pre-filtering the 2-d points in field image that are far from the epipolar line corresponding to the TargetPoint
80
81        % Computing epipolar lines
82        EpipolarLinePara = F*[ TargetPointPix(:,IndexTargetPointPix); 1];
83        NormalizedEpipolarLinePara = EpipolarLinePara / norm(EpipolarLinePara(1:2));   
84
85        % Find the RePojection of the Target Point in Field-Camera
86        TargetPoint3D = ModelInfoTarget.Model.PlaneParaInfo.Position3DFited(:,INDTarget(:,IndexTargetPointPix));    % Size 3x1
87        PointDepth(IndexTargetPointPix) = ModelInfoTarget.Model.PlaneParaInfo.FitDepth(INDTarget(:,IndexTargetPointPix));
88        POriReproj =  defaultPara.InrinsicK2*(Pair.R*( TargetPoint3D) + Pair.T);
89        POriReproj = POriReproj(1:2,:)./POriReproj(3,:);
90
91        % Modify the POriReproj to the closest point in the image
92        % Since ReProjection of the Target Point might not be in size the Field image range
93        POriReproj = Point2ImageRange(POriReproj, NormalizedEpipolarLinePara, ScaleImg); % some error/ Min fixed later//Solved 6/28
94
95        % Finding Field points with distances to the epipolar line smaller than EpipolarThre*max(ScaleImg)
96        Mask = abs(NormalizedEpipolarLinePara' *  [ FieldPointPix; ones(1, size(FieldPointPix,2))]) <= EpipolarThre*max(ScaleImg);
97        INDFieldCandidate = INDField(Mask);
98        FieldPoint3D = AllFieldPoint3D(:,INDFieldCandidate);    % 3 x NPoints % inFeild coordinate
99   
100   
101    % 2) Pre-filtering the 3-d points in field image that are far from the ray connecting target point to target image-origin.
102   
103        % Find the distance of all the points from the target-origi ray
104        if ~isempty(FieldPoint3D)
105                u = TargetPoint3D' * FieldPoint3D / (TargetPoint3D' * TargetPoint3D);
106                distances = sqrt( sum( ( repmat(u,[3,1]) .* repmat(TargetPoint3D,[1,size(FieldPoint3D,2)]) - FieldPoint3D).^2 ,1));
107                occludingPointsMask = distances < MaxFieldPointDist2Ray;
108                distances = distances(occludingPointsMask);
109                [distances SortIND ] = sort(distances);
110                % update the new Field Candidate Points
111                FieldPoint3D = FieldPoint3D(:,occludingPointsMask);
112                FieldPoint3D = FieldPoint3D(:,SortIND);
113                INDFieldCandidate = INDFieldCandidate(occludingPointsMask);
114                INDFieldCandidate = INDFieldCandidate(SortIND);
115        end   
116   
117    % ============ Min's New trial =================
118    % 3) Check triangals using FieldPoint3D are occluing the Ray or Not
119        [ SubV SubH] = ind2sub(DepthScale, INDFieldCandidate);
120%     length(INDFieldCandidate)
121    i = 1;
122    breakFlag = 0;
123        while i <=length(INDFieldCandidate)     && breakFlag == 0
124%         i
125                Ray =  TargetPoint3D./norm(TargetPoint3D);
126                if (~(SubV(i) == 1 || SubH(i) == 1)     && breakFlag ==0)
127                % check top left square
128                        if ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i)-1, SubH(i)-1) == ...
129                                ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i), SubH(i))
130                                % upper right
131                                FaceSetID1 = INDFieldCandidate(i) - DepthScale(1) - 1;
132                                FaceSetID2 = INDFieldCandidate(i) - 1;
133                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)-1 SubV(i)-1],[SubH(i) SubH(i)-1 SubH(i)]);
134                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
135                                        Ray);
136                                OccluStorageScript; % min added script
137                                % lower left
138                                FaceSetID1 = INDFieldCandidate(i) - DepthScale(1) - 1;
139                                FaceSetID2 = INDFieldCandidate(i) - DepthScale(1);
140                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)-1 SubV(i)],[SubH(i) SubH(i)-1 SubH(i)-1]);
141                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
142                                        Ray);
143                                OccluStorageScript; % min added script
144                        else
145                                % lower right
146                                FaceSetID1 = INDFieldCandidate(i) - DepthScale(1);
147                                FaceSetID2 = INDFieldCandidate(i) - 1;
148                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i) SubV(i)-1],[SubH(i) SubH(i)-1 SubH(i)]);
149                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
150                                        Ray);
151                                OccluStorageScript; % min added script
152                        end
153                end
154
155                if (~(SubV(i) == 1 || SubH(i) == DepthScale(2)) && breakFlag ==0)
156                % check top right square
157                        if ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i)-1, SubH(i)) == ...
158                                ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i), SubH(i)+1)
159                                % lower left
160                                FaceSetID1 = INDFieldCandidate(i) - 1;
161                                FaceSetID2 = INDFieldCandidate(i) + DepthScale(1);
162                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)-1 SubV(i)],[SubH(i) SubH(i) SubH(i)+1]);
163                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
164                                        Ray);
165                                OccluStorageScript; % min added script
166                        else
167                                % upper left
168                                FaceSetID1 = INDFieldCandidate(i) - 1;
169                                FaceSetID2 = INDFieldCandidate(i) + DepthScale(1) - 1;
170                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)-1 SubV(i)-1],[SubH(i) SubH(i) SubH(i)+1]);
171                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
172                                        Ray);
173                                OccluStorageScript; % min added script
174                                % lower right
175                                FaceSetID1 = INDFieldCandidate(i) + DepthScale(1) - 1;
176                                FaceSetID2 = INDFieldCandidate(i) + DepthScale(1);
177                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)-1 SubV(i)],[SubH(i) SubH(i)+1 SubH(i)+1]);
178                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
179                                        Ray);
180                                OccluStorageScript; % min added script
181
182                        end
183                end
184               
185                if (~(SubV(i) == DepthScale(1) || SubH(i) == 1) && breakFlag ==0)
186                % check bottom left square
187                        if ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i), SubH(i)-1) == ...
188                                ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i)+1, SubH(i))
189                                % upper right
190                                FaceSetID1 = INDFieldCandidate(i) - DepthScale(1);
191                                FaceSetID2 = INDFieldCandidate(i) + 1;
192                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i) SubV(i)+1],[SubH(i) SubH(i)-1 SubH(i)]);
193                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
194                                        Ray);
195                                OccluStorageScript; % min added script
196
197                        else
198                                % upper left
199                                FaceSetID1 = INDFieldCandidate(i) - DepthScale(1);
200                                FaceSetID2 = INDFieldCandidate(i) - DepthScale(1) + 1;
201                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i) SubV(i)+1],[SubH(i) SubH(i)-1 SubH(i)-1]);
202                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
203                                        Ray);
204                                OccluStorageScript; % min added script
205                                % lower right
206                                FaceSetID2 = INDFieldCandidate(i) - DepthScale(1) + 1;
207                                FaceSetID1 = INDFieldCandidate(i) + 1;
208                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)+1 SubV(i)+1],[SubH(i) SubH(i)-1 SubH(i)]);
209                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
210                                        Ray);
211                                OccluStorageScript; % min added script
212                        end
213                end
214
215                if (~(SubV(i) == DepthScale(1) || SubH(i) == DepthScale(2)) && breakFlag ==0)
216                % check bottom left square
217                        if ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i), SubH(i)) == ...
218                                ModelInfoField.Model.PlaneParaInfo.SupOri(SubV(i)+1, SubH(i)+1)
219                                % upper right
220                                FaceSetID2 = INDFieldCandidate(i) + DepthScale(1) + 1;
221                                FaceSetID1 = INDFieldCandidate(i) + DepthScale(1);
222                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)+1 SubV(i)],[SubH(i) SubH(i)+1 SubH(i)+1]);
223                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
224                                        Ray);
225                                OccluStorageScript; % min added script
226                                % lower left
227                                FaceSetID2 = INDFieldCandidate(i) + 1;
228                                FaceSetID1 = INDFieldCandidate(i) + DepthScale(1) + 1;
229                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)+1 SubV(i)+1],[SubH(i) SubH(i) SubH(i)+1]);
230                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
231                                        Ray);
232                                OccluStorageScript; % min added script
233                        else
234                                % upper left
235                                FaceSetID2 = INDFieldCandidate(i) + 1;
236                                FaceSetID1 = INDFieldCandidate(i) + DepthScale(1);
237                                TriInd = sub2ind(DepthScale,[ SubV(i) SubV(i)+1 SubV(i)],[SubH(i) SubH(i) SubH(i)+1]);
238                                [ CombinePara OccluPointDepth] = SolveOccluPoint( AllFieldPoint3D(:,TriInd),...
239                                        Ray);
240                                OccluStorageScript; % min added script
241                        end
242        end
243%     i
244    i = i + 1;
245    end % end of second for loop
246        % =====================================
247% IndexTargetPointPix
248end % end of first for loop
249disp([ '        ' num2str( toc(FindOccluPointTime)) ' seconds']);
250return;
Note: See TracBrowser for help on using the repository browser.