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 | % */
|
---|
39 | function [SupMerg MedSupMerg] = FitPlaneLaser_CoPlane(... |
---|
40 | Sup,MedSup,depthMap,Ray,MedRay,maskSky,maskG,Algo,k, CornerList, OccluList,previoslyStored); |
---|
41 | % This function runs the RMF for the plane parameter of each superpixels |
---|
42 | % output: |
---|
43 | % newSup : newSuperpixel index |
---|
44 | % PlanePara : Plane Parameter corresponding to the new Superpuixel index |
---|
45 | |
---|
46 | displayFlag = false; |
---|
47 | |
---|
48 | global GeneralDataFolder ScratchDataFolder LocalFolder ClusterExecutionDirectory... |
---|
49 | ImgFolder VertYNuPatch VertYNuDepth HoriXNuPatch HoriXNuDepth a_default b_default Ox_default Oy_default... |
---|
50 | Horizon_default filename batchSize NuRow_default SegVertYSize SegHoriXSize WeiBatchSize PopUpVertY PopUpHoriX taskName; |
---|
51 | |
---|
52 | TrialName ='LaserL1PosiStickAlign' |
---|
53 | Mult =1; |
---|
54 | vari =0.1; |
---|
55 | shift = 10; |
---|
56 | StickHori = 1; |
---|
57 | StickVert = 10; |
---|
58 | HoriConf = 0; |
---|
59 | VertConf = 0; |
---|
60 | mapVert = logspace(VertConf,0,VertYNuDepth) % modeling the gravity prior |
---|
61 | mapHori = [logspace(HoriConf,0,round(HoriXNuDepth/2)) fliplr(logspace(HoriConf,0,HoriXNuDepth-round(HoriXNuDepth/2)))] % modeling the gravity prior |
---|
62 | gravity =true; |
---|
63 | threhConer = 0.5; |
---|
64 | threhJump = 1.2; |
---|
65 | ceiling = 0*VertYNuDepth; |
---|
66 | % get rid of the error point and sky point in the depthMap |
---|
67 | % set every depth bigger than 80meter as error point and sky point |
---|
68 | %CleanedDepthMap = (depthMapif ~previoslyStored >80).*medfilt2(depthMap,[4 4])+(depthMap<=80).*depthMap; |
---|
69 | CleanedDepthMap = depthMap; |
---|
70 | %if ~learned |
---|
71 | % CleanedDepthMap(depthMap>80) = NaN; % don't clean the point >80 sometimes it occlusion |
---|
72 | %end |
---|
73 | Posi3D = im_cr2w_cr(CleanedDepthMap,permute(Ray,[2 3 1])); |
---|
74 | |
---|
75 | %previoslyStored = true; |
---|
76 | if ~previoslyStored |
---|
77 | |
---|
78 | NewMap = [rand(max(Sup(:)),3); [0 0 0]]; |
---|
79 | % Clean the Sup |
---|
80 | %maskSky = imdilate(maskSky, strel('disk', 3) ); |
---|
81 | %[Sup,MedSup]=CleanedSup(Sup,MedSup,maskSky); |
---|
82 | maskSky = Sup == 0; |
---|
83 | maskSkyEroded = imerode(maskSky, strel('disk', 4) ); |
---|
84 | SupEpand = ExpandSup2Sky(Sup,maskSkyEroded); |
---|
85 | NuPatch = HoriXNuDepth*VertYNuDepth-sum(maskSky(:)); |
---|
86 | |
---|
87 | |
---|
88 | % find out the neighbors |
---|
89 | %[nList]=GenSupNeighbor(Sup,maskSky); % Cell array : [ Supindex, Neighbor List of SupIndex] |
---|
90 | NuSup = setdiff(unique(Sup)',0); |
---|
91 | NuSup = sort(NuSup); |
---|
92 | NuSupSize = size(NuSup,2); |
---|
93 | |
---|
94 | % Sup index and planeParameter index inverse map |
---|
95 | Sup2Para = sparse(1,max(Sup(:))); |
---|
96 | Sup2Para(NuSup) = 1:NuSupSize; |
---|
97 | |
---|
98 | % Generate the Matrix for MRF |
---|
99 | tic |
---|
100 | PosiM = sparse(0,0); |
---|
101 | RayMd = sparse(0,0); |
---|
102 | RayAllM = sparse(0,0); |
---|
103 | RayMtilt = sparse(0,0); |
---|
104 | RayMCent = sparse(0,0); |
---|
105 | DepthInverseMCent = []; |
---|
106 | DepthInverseM = []; |
---|
107 | YPointer = []; |
---|
108 | beta = []; |
---|
109 | EmptyIndex = []; |
---|
110 | for i = NuSup |
---|
111 | mask = Sup ==i; |
---|
112 | RayAllM = blkdiag( RayAllM, Ray(:,mask)'); |
---|
113 | [yt x] = find(mask); |
---|
114 | CenterX = round(median(x)); |
---|
115 | CenterY = round(median(yt)); |
---|
116 | YPointer = [YPointer; CenterY >= ceiling]; |
---|
117 | mask(isnan(CleanedDepthMap)) = false; |
---|
118 | SupNuPatch(i) = sum(mask(:)); |
---|
119 | % if sum(mask(:)) < 5 |
---|
120 | % EmptyIndex = [EmptyIndex; i]; |
---|
121 | % mask(:) = false; |
---|
122 | % end |
---|
123 | % find center point |
---|
124 | [yt x] = find(mask); |
---|
125 | CenterX = round(median(x)); |
---|
126 | CenterY = round(median(yt)); |
---|
127 | |
---|
128 | if ~all(mask(:)==0) |
---|
129 | if gravity |
---|
130 | Pa2CenterRatio = median(CleanedDepthMap(mask))./CleanedDepthMap(mask); |
---|
131 | if sum(mask(:)) > 0 |
---|
132 | RayMtilt = blkdiag(RayMtilt, ... |
---|
133 | ( Pa2CenterRatio(:,[1 1 1]).*repmat(Ray(:,CenterY,CenterX)',[ SupNuPatch(i) 1])- Ray(:,mask)')); |
---|
134 | else |
---|
135 | RayMtilt = blkdiag(RayMtilt, Ray(:,mask)'); |
---|
136 | end |
---|
137 | RayMCent = blkdiag(RayMCent, Ray(:,CenterY,CenterX)'*SupNuPatch(i)*mapVert(CenterY)*mapHori(CenterX)); |
---|
138 | PosiM = blkdiag(PosiM,Posi3D(:,mask)'); |
---|
139 | RayMd = blkdiag(RayMd,Ray(:,mask)'.*repmat( mapVert(yt)',[1 3]).*repmat( mapHori(x)',[1 3])); |
---|
140 | DepthInverseMCent = [DepthInverseMCent; 1./median(CleanedDepthMap(mask)).*SupNuPatch(i)* mapVert(CenterY)'*mapHori(CenterX)]; |
---|
141 | DepthInverseM = [DepthInverseM; 1./CleanedDepthMap(mask).* mapVert(yt)'.*mapHori(x)']; |
---|
142 | else |
---|
143 | Pa2CenterRatio = CleanedDepthMap(CenterY,CenterX)./CleanedDepthMap(mask); |
---|
144 | if sum(mask(:)) > 0 |
---|
145 | RayMtilt = blkdiag(RayMtilt, ... |
---|
146 | ( Pa2CenterRatio(:,[1 1 1]).*repmat(Ray(:,CenterY,CenterX)',[ SupNuPatch(i) 1])- Ray(:,mask)')); |
---|
147 | else |
---|
148 | RayMtilt = blkdiag(RayMtilt, Ray(:,mask)'); |
---|
149 | end |
---|
150 | RayMCent = blkdiag(RayMCent, Ray(:,CenterY,CenterX)'*SupNuPatch(i)); |
---|
151 | PosiM = blkdiag(PosiM,Posi3D(:,mask)'); |
---|
152 | RayMd = blkdiag(RayMd,Ray(:,mask)'); |
---|
153 | DepthInverseMCent = [DepthInverseMCent; 1./median(CleanedDepthMap(mask)).*SupNuPatch(i)]; |
---|
154 | DepthInverseM = [DepthInverseM; 1./CleanedDepthMap(mask)]; |
---|
155 | end |
---|
156 | else |
---|
157 | RayMtilt = blkdiag(RayMtilt, Ray(:,mask)'); |
---|
158 | RayMCent = blkdiag(RayMCent, Ray(:,mask)'); |
---|
159 | PosiM = blkdiag(PosiM, Posi3D(:,mask)'); |
---|
160 | RayMd = blkdiag(RayMd, Ray(:,mask)'); |
---|
161 | end |
---|
162 | end |
---|
163 | YPointer(YPointer==0) = -1; |
---|
164 | % buliding CoPlane Matrix========================================================================= |
---|
165 | BounaryPHori = conv2(Sup,[1 -1],'same') ~=0; |
---|
166 | BounaryPHori(:,end) = 0; |
---|
167 | BounaryPVert = conv2(Sup,[1; -1],'same') ~=0; |
---|
168 | BounaryPVert(end,:) = 0; |
---|
169 | ClosestNList = [ Sup(find(BounaryPHori==1)) Sup(find(BounaryPHori==1)+VertYNuDepth);... |
---|
170 | Sup(find(BounaryPVert==1)) Sup(find(BounaryPVert==1)+1)]; |
---|
171 | ClosestNList = sort(ClosestNList,2); |
---|
172 | ClosestNList = unique(ClosestNList,'rows'); |
---|
173 | ClosestNList(ClosestNList(:,1) == 0,:) = []; |
---|
174 | BoundaryAll = BounaryPHori + BounaryPHori(:,[1 1:(end-1)])... |
---|
175 | +BounaryPVert + BounaryPVert([1 1:(end-1)],:); |
---|
176 | BoundaryAll([1 end],:) = 1; |
---|
177 | BoundaryAll(:,[1 end]) = 1; |
---|
178 | NuNei = size(ClosestNList,1); |
---|
179 | CoPM1 = sparse(0,3*NuSupSize); |
---|
180 | CoPM2 = sparse(0,3*NuSupSize); |
---|
181 | WeiCoP = []; |
---|
182 | for i = 1: NuNei |
---|
183 | % if ~CornerList(i) |
---|
184 | mask = Sup == ClosestNList(i,1); |
---|
185 | SizeMaskAll = sum(mask(:)); |
---|
186 | [y x] = find(mask); |
---|
187 | CenterX = round(median(x)); |
---|
188 | CenterY = round(median(y)); |
---|
189 | y = find(mask(:,CenterX)); |
---|
190 | if ~isempty(y) |
---|
191 | CenterY = round(median(y)); |
---|
192 | end |
---|
193 | % [yt x] = find(mask); |
---|
194 | % CenterX = round(median(x)); |
---|
195 | % CenterY = round(median(yt)); |
---|
196 | |
---|
197 | temp1 = sparse(1, 3*NuSupSize); |
---|
198 | temp2 = sparse(1, 3*NuSupSize); |
---|
199 | temp1(:,(Sup2Para(ClosestNList(i,1))*3-2): Sup2Para(ClosestNList(i,1))*3) = Ray(:,CenterY,CenterX)'; |
---|
200 | temp2(:,(Sup2Para(ClosestNList(i,2))*3-2): Sup2Para(ClosestNList(i,2))*3) = Ray(:,CenterY,CenterX)'; |
---|
201 | CoPM1 = [CoPM1; temp1]; |
---|
202 | CoPM2 = [CoPM2; temp2]; |
---|
203 | tempWeiCoP = [SizeMaskAll]; |
---|
204 | |
---|
205 | mask = Sup == ClosestNList(i,2); |
---|
206 | SizeMaskAll = sum(mask(:)); |
---|
207 | [y x] = find(mask); |
---|
208 | CenterX = round(median(x)); |
---|
209 | CenterY = round(median(y)); |
---|
210 | y = find(mask(:,CenterX)); |
---|
211 | if ~isempty(y) |
---|
212 | CenterY = round(median(y)); |
---|
213 | end |
---|
214 | % [yt x] = find(mask); |
---|
215 | % CenterX = round(median(x)); |
---|
216 | % CenterY = round(median(yt)); |
---|
217 | |
---|
218 | temp1 = sparse(1, 3*NuSupSize); |
---|
219 | temp2 = sparse(1, 3*NuSupSize); |
---|
220 | temp1(:,(Sup2Para(ClosestNList(i,1))*3-2): Sup2Para(ClosestNList(i,1))*3) = Ray(:,CenterY,CenterX)'; |
---|
221 | temp2(:,(Sup2Para(ClosestNList(i,2))*3-2): Sup2Para(ClosestNList(i,2))*3) = Ray(:,CenterY,CenterX)'; |
---|
222 | CoPM1 = [CoPM1; temp1]; |
---|
223 | CoPM2 = [CoPM2; temp2]; |
---|
224 | tempWeiCoP = [tempWeiCoP; SizeMaskAll]; |
---|
225 | WeiCoP = [WeiCoP; tempWeiCoP]; |
---|
226 | % end |
---|
227 | end%========================================================================================================= |
---|
228 | |
---|
229 | % find the boundary point that might need to be stick ot each other========================================== |
---|
230 | HoriStickM = sparse(1,2+3*NuSupSize); |
---|
231 | for i = find(BounaryPHori==1)' |
---|
232 | j = i+VertYNuDepth; |
---|
233 | if Sup(i) == 0 || Sup(j) == 0 |
---|
234 | continue; |
---|
235 | end |
---|
236 | % if ~OccluList(sum( ClosestNList == repmat(sort([Sup(i) Sup(j)]), [NuNei 1]),2) == 2) |
---|
237 | Target(1) = Sup2Para(Sup(i)); |
---|
238 | Target(2) = Sup2Para(Sup(j)); |
---|
239 | rayBoundary(:,1) = Ray(:,i); |
---|
240 | rayBoundary(:,2) = Ray(:,i); |
---|
241 | betaTemp = StickHori;%*(DistStickLengthNormWei.^2)*beta(Target(I)); |
---|
242 | temp = sparse(3,NuSupSize); |
---|
243 | temp(:,Target(1)) = rayBoundary(:,1); |
---|
244 | temp(:,Target(2)) = -rayBoundary(:,2); |
---|
245 | HoriStickM = [HoriStickM; [sort([Sup(i) Sup(j)]) betaTemp*temp(:)']]; |
---|
246 | % end |
---|
247 | end |
---|
248 | VertStickM = sparse(1,2+3*NuSupSize); |
---|
249 | for i = find(BounaryPVert==1)' |
---|
250 | j = i+1; |
---|
251 | if Sup(i) == 0 || Sup(j) == 0 |
---|
252 | continue; |
---|
253 | end |
---|
254 | % if ~OccluList(sum( ClosestNList == repmat(sort([Sup(i) Sup(j)]), [NuNei 1]),2) == 2) |
---|
255 | Target(1) = Sup2Para(Sup(i)); |
---|
256 | Target(2) = Sup2Para(Sup(j)); |
---|
257 | rayBoundary(:,1) = Ray(:,i); |
---|
258 | rayBoundary(:,2) = Ray(:,i); |
---|
259 | betaTemp = StickVert;%*(DistStickLengthNormWei.^2)*beta(Target(I)); |
---|
260 | temp = sparse(3,NuSupSize); |
---|
261 | temp(:,Target(1)) = rayBoundary(:,1); |
---|
262 | temp(:,Target(2)) = -rayBoundary(:,2); |
---|
263 | VertStickM = [VertStickM; [sort([Sup(i) Sup(j)]) betaTemp*temp(:)']]; |
---|
264 | % end |
---|
265 | end |
---|
266 | |
---|
267 | % 2nd: Associate Plane Parameter term and Interactive Co-Planar term |
---|
268 | opt = sdpsettings('solver','sedumi'); |
---|
269 | ParaPPCP = sdpvar(3*NuSupSize,1); |
---|
270 | F = set(RayAllM*ParaPPCP >=0.01); |
---|
271 | %F = set(ParaPPCP(3*(1:NuSupSize)-1)<=0) + set(RayAllM*ParaPPCP >=0.01); |
---|
272 | % First fit the plane to find the estimated plane parameters |
---|
273 | % If PosiM contain NaN data the corresponding Plane Parameter will also be NaN |
---|
274 | %solvesdp([], norm( spdiags(ConfRayMean',0,size(RayMeanM,1),size(RayMeanM,1))*(RayMeanM*ParaOri(:) - RayMeanM*ParaPPCP),1)... |
---|
275 | %solvesdp([],norm( PosiM*ParaPPCP-ones(size(PosiM,1),1),1)... |
---|
276 | %solvesdp(F,norm( RayMCent*ParaPPCP - DepthInverseMCent,1) ... |
---|
277 | solvesdp(F,norm( RayMd*ParaPPCP - DepthInverseM,1)... |
---|
278 | +norm((CoPM1 - CoPM2)*ParaPPCP,1), opt); |
---|
279 | % +norm(VertStickM(:,3:end)*ParaPPCP,1)+ norm(HoriStickM(:,3:end)*ParaPPCP,1), opt); |
---|
280 | % +0.01*norm( RayMtilt*ParaPPCP,1)... |
---|
281 | % +norm(spdiags(WeiCoP,0,size(CoPM1,1),size(CoPM1,1))*(CoPM1 - CoPM2)*ParaPPCP,1)... |
---|
282 | % +norm( (CoPM1 - CoPM2)*ParaPPCP,1 )... |
---|
283 | ParaPPCP = double(ParaPPCP); |
---|
284 | ParaPPCP = reshape(ParaPPCP,3,[]); |
---|
285 | any(any(isnan(ParaPPCP))) |
---|
286 | |
---|
287 | % porject the ray on planes to generate the ProjDepth |
---|
288 | FitDepthPPCP = 80*ones(1,VertYNuDepth*HoriXNuDepth); |
---|
289 | FitDepthPPCP(~maskSkyEroded) = (1./sum(ParaPPCP(:,Sup2Para(SupEpand(~maskSkyEroded))).*Ray(:,~maskSkyEroded),1))'; |
---|
290 | FitDepthPPCP = reshape(FitDepthPPCP,VertYNuDepth,[]); |
---|
291 | [Position3DFitedPPCP] = im_cr2w_cr(FitDepthPPCP,permute(Ray,[2 3 1])); |
---|
292 | [VrmlName] = vrml_test_faceset_goodSkyBoundary( filename{k}, Position3DFitedPPCP, FitDepthPPCP, permute(Ray,[2 3 1]), 'NoBLaser', ... |
---|
293 | [], [], 0, maskSkyEroded, maskG, 1, 0, a_default, b_default, Ox_default, Oy_default); |
---|
294 | system(['gzip -9 -c ' ScratchDataFolder '/vrml/' VrmlName ' > ' ScratchDataFolder '/vrml/' VrmlName '.gz']); |
---|
295 | delete([ScratchDataFolder '/vrml/' VrmlName]); |
---|
296 | %return; |
---|
297 | %======================================================================================================================================= |
---|
298 | % Calculate Occlusion metrics |
---|
299 | [SpatialJump] = DecideOcclusion(VertStickM,HoriStickM,ParaPPCP,ClosestNList);%[VertStickM(:,1:2) VertStickM(:,3:end)*ParaPPCP; HoriStickM(:,1:2) HoriStickM(:,3:end)*ParaPPCP]; |
---|
300 | OccluList = SpatialJump > log(threhJump); |
---|
301 | [SpatialDisMask] = ShowCorner(MedSup,[ClosestNList OccluList]); |
---|
302 | |
---|
303 | % CalCulate Corner Metrics |
---|
304 | NormParaPPCP = ParaPPCP./repmat(norms(ParaPPCP),[3 1]); |
---|
305 | ClosestNListCorner = [ClosestNList sum(abs(NormParaPPCP(:,Sup2Para(ClosestNList(:,1))) - NormParaPPCP(:,Sup2Para(ClosestNList(:,2)))),1)']; |
---|
306 | CornerList = (ClosestNListCorner(:,3)/3)>threhConer ; |
---|
307 | [CornerMask] = ShowCorner(MedSup,[ClosestNListCorner(:,1:2) CornerList ]); |
---|
308 | save([ScratchDataFolder '/data/temp/List' num2str(k) '.mat'],'ClosestNList','ClosestNListCorner','CornerMask','SpatialJump','SpatialDisMask','OccluList','CornerList'); |
---|
309 | |
---|
310 | end |
---|
311 | return; |
---|