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 ReportPlaneParaMRF_Conditioned(Default, Sup, SupOri,depthMap,VarMapRaw, ... |
---|
40 | RayOri, Ray, SupNeighborTable, maskSky,maskG, MultiScaleSupTable, ... |
---|
41 | StraightLineTable); |
---|
42 | % This function runs the RMF over the plane parameter of each superpixels |
---|
43 | |
---|
44 | % Input: |
---|
45 | % Default- all default parameters |
---|
46 | % Sup - Superpixel, in 2-d matrix. |
---|
47 | % SupOri - original superpixel before cleaning |
---|
48 | % DepthMap - DepthMap from Learning or coarse measuring |
---|
49 | % VarMap - confidance for DepthMap (latest algo VarMap used simply one since the learning not generalized well) |
---|
50 | % RayOri - Rays without stiching to the boundary. |
---|
51 | % Ray - rays after stiching |
---|
52 | % SupNeighborTable - A lookup table for superpixels' neighbors (2-d sparse matrix) |
---|
53 | % maskSky - Sky mask |
---|
54 | % maskG - Ground mask |
---|
55 | % MultiScaleSupTable - multiple scale segmentation, used to define the weights betwenn the stiching terms. |
---|
56 | % StraightLineTable - straight line stiching, (not used in this version, but might be very usefulif implemented). |
---|
57 | |
---|
58 | % ------------------------------------------Finish Input Definition ------------------------------------------------ |
---|
59 | % Parameter setting |
---|
60 | % 1) Functional parameters |
---|
61 | FlagRemoveVerticalSupport = 1; % Enable Removing Vertical Support Depth at the second inference (risky if occluded) |
---|
62 | FlagEnableVarMap = 0; % Disabled |
---|
63 | gravity =true; % if true, apply the HoriConf and VertConf linear scale weight |
---|
64 | CoPST = true; % if true, apply the Straight line prior as the Co-Planar constrain |
---|
65 | % ============= Magic Number ============= |
---|
66 | StickHori = 5; %0.1; % sticking power in horizontal direction |
---|
67 | StickVert = 10; % sticking power in vertical direction |
---|
68 | Center = 2; % Co-Planar weight at the Center of each superpixel |
---|
69 | HoriConf = 1; % set the confidant of the learned depth at the middle in Horizontal direction of the image |
---|
70 | VertConf = 0.01; % set the confidant of the learned depth at the top of the image |
---|
71 | BandWith = 1; % Nov29 1 Nov30 0.1 check result change to 0.1 12/1 0.1 lost detail |
---|
72 | ClosestDist = 1; % define the closest distance that the MRF allows |
---|
73 | FarestDist = 1.5*median(depthMap(:)); % tried on university % nogood effect but keep it since usually it the rangelike this % change to 1.5 for church |
---|
74 | % The hand-tuned 14 dimensional vector correspose to the weights of 14 multiple segmentation |
---|
75 | if ~isempty(MultiScaleSupTable) |
---|
76 | MultiScaleFlag = true; % multiple segmentaition hypothesis |
---|
77 | WeiV = 2*ones(1,size(MultiScaleSupTable,2)-1); |
---|
78 | else |
---|
79 | MultiScaleFlag = false; |
---|
80 | WeiV = 1; |
---|
81 | end |
---|
82 | WeiV(1,1:2:end) = 6; % emphasize the middle scale three times smaller than large scale |
---|
83 | WeiV = WeiV./sum(WeiV);% normalize if pair of superpixels have same index in all the scale, their weight will be 10 |
---|
84 | ShiftStick = -.1; % between -1 and 0, more means more smoothing. |
---|
85 | ShiftCoP = -.5; % between -1 and 0, more means more smoothing. |
---|
86 | % ========================================================================= |
---|
87 | % If you go up in the vertical direction in the image, the weights change in the vertical direction. |
---|
88 | groundThreshold = cos([ zeros(1, Default.VertYNuDepth - ceil(Default.VertYNuDepth/2)+10) ... |
---|
89 | linspace(0,15,ceil(Default.VertYNuDepth/2)-10)]*pi/180); |
---|
90 | % v1 15 v2 20 too big v3 20 to ensure non misclassified as ground. |
---|
91 | % verticalThreshold = cos(linspace(5,55,Default.VertYNuDepth)*pi/180); % give a vector of size 55 in top to down : |
---|
92 | verticalThreshold = cos([ 5*ones(1,Default.VertYNuDepth - ceil(Default.VertYNuDepth/2)) ... |
---|
93 | linspace(5,55,ceil(Default.VertYNuDepth/2))]*pi/180); |
---|
94 | % give a vector of size 55 in top to down : |
---|
95 | % 50 means suface norm away from y axis more than 50 degree |
---|
96 | % ========================================================================= |
---|
97 | ceiling = 0*Default.VertYNuDepth; % set the position of the ceiling, related to No plane coming back constrain % changed for newchurch |
---|
98 | % ^^^^^This number means the lowest row that planes may go back just like a ceiling |
---|
99 | % ============= End of Magic Number ====== |
---|
100 | |
---|
101 | % 2) Opimization parameters |
---|
102 | % ============== parameters for the decomposition problem |
---|
103 | % Optimal parameters for current code: |
---|
104 | % For one machine: Sedumi: (1,1) Lpsolve:(3,1) |
---|
105 | % Multiple machines: Sedumi: (4,1) LPsolve:(3,1) |
---|
106 | % lpsolve running time is 22 seconds for (4,2) arrangement; but numerical |
---|
107 | % accuracy needs to be resolved first. |
---|
108 | XNuDecompose = 1; % up to 3 is stable |
---|
109 | YNuDecompose = 1; |
---|
110 | % ============ parameters for the decomposition problem |
---|
111 | solverVerboseLevel = 0; % set to 1 if you need msg out when solving optimization problem |
---|
112 | NoSecondStep = 0; % set to 1 if only want first level opt = no second step of vertical and horizontal objective |
---|
113 | |
---|
114 | % 3) Debug and evaluation paramters |
---|
115 | ExtractRelationInfo = 1; % set to 1 if need to storage the coplanar and connectivity weight to analyze |
---|
116 | |
---|
117 | % 4) Rendering parameters |
---|
118 | GridFlag = 0; % set to 1 if the wrl need grid of the triangle overlay |
---|
119 | |
---|
120 | % ==========================End of parameter setting ==================================== |
---|
121 | inferenceTime = tic; |
---|
122 | fprintf(['\n : Building Matrices.... ']); |
---|
123 | |
---|
124 | % ======= intermediant Data ===================================================================== |
---|
125 | % confidence of the supporting depth changes according to the row and column ========= |
---|
126 | mapVert = linspace(VertConf,1,Default.VertYNuDepth); % modeling the gravity prior: the topper the row is the lower the confidence of supporting depths |
---|
127 | mapHori = [linspace(HoriConf,1,round(Default.HoriXNuDepth/2)) fliplr(linspace(HoriConf,1,Default.HoriXNuDepth-round(Default.HoriXNuDepth/2)))]; |
---|
128 | % the more peripheral column the lower the confidence of supporting depths |
---|
129 | % ==================================================================================== |
---|
130 | % assign the confidance of the supporting depths |
---|
131 | if FlagEnableVarMap |
---|
132 | VarMap = zeros( size(VarMapRaw)); |
---|
133 | else |
---|
134 | VarMap = VarMapRaw; |
---|
135 | end |
---|
136 | CleanedDepthMap = depthMap; |
---|
137 | CleanedDepthMap(depthMap>FarestDist) = NaN; % set the supporting depths ti NaN means that depth is not effective in OPT, so we simply ignore the depth which is too far |
---|
138 | Posi3D = im_cr2w_cr(CleanedDepthMap,permute(RayOri,[2 3 1])); % given the depth and ray as input, calculate the 3-d coordinate at each point. |
---|
139 | |
---|
140 | % Clean the Sup near sky |
---|
141 | maskSky = Sup == 0; |
---|
142 | maskSkyEroded = imerode(maskSky, strel('disk', 4) ); |
---|
143 | SupEpand = ExpandSup2Sky(Sup,maskSkyEroded); % extend the sky using the cloesest Sup index, which mean the some part of the sky will be include in the sup but their supporting depth will not be used. |
---|
144 | |
---|
145 | NuSup = setdiff(unique(Sup)',0); % unique index of sup (not including sky index which is '0') |
---|
146 | NuSupSize = size(NuSup,2); % number of uniquew sup |
---|
147 | |
---|
148 | % Sup index and planeParameter index inverse map ======= useful tool====== |
---|
149 | % since sup index is not continuous, but the Parameter index over sup is continuous. Sup2Para is the maping stand for sup index to Parameter index |
---|
150 | Sup2Para = sparse(1,max(Sup(:))); |
---|
151 | Sup2Para(NuSup) = 1:NuSupSize; |
---|
152 | % ======================================================================== |
---|
153 | |
---|
154 | % =================please ignore here since StraightLineTable = [], NuLine = 0 ======== |
---|
155 | % constructiion of the Straight line prior matrix Will be add in the CoPlane matrix |
---|
156 | NuLine = size(StraightLineTable,2); |
---|
157 | CoPSTList = []; |
---|
158 | |
---|
159 | % effectively not running the straight line constraint here. |
---|
160 | for i = 1:NuLine |
---|
161 | L = StraightLineTable{i}; |
---|
162 | X = L(1:(end-1))'; |
---|
163 | Y = L(2:end)'; |
---|
164 | if isempty(X) |
---|
165 | continue; |
---|
166 | end |
---|
167 | for j = 1:size(X,1) |
---|
168 | if X(j)~=Y(j) |
---|
169 | CoPSTList = [CoPSTList; X(j) Y(j)]; |
---|
170 | end |
---|
171 | end |
---|
172 | end |
---|
173 | % Please ignore =============================================================== |
---|
174 | |
---|
175 | % ============end of generating intermediant Data =============================================== |
---|
176 | |
---|
177 | % Generate the Matrix for MRF ============================================================== |
---|
178 | PosiM = sparse(0,0); % supporting depth times ray terms |
---|
179 | VarM = sparse(0,0); % confidence terms |
---|
180 | RayAllM = sparse(0,0); % collection of every ray |
---|
181 | |
---|
182 | % keep record of the center of the sup is lower then the ceiling |
---|
183 | YPointer = []; % set to one if lower then the ceiling |
---|
184 | YPosition = []; % keep a record of the row of the center of the sup, useful when enforcing vertical and horizontal constrain |
---|
185 | |
---|
186 | for i = NuSup |
---|
187 | mask = SupEpand ==i; % include the Ray that will be use to expand the NonSky |
---|
188 | RayAllM = blkdiag( RayAllM, Ray(:,mask)'); % RayAllM will be putting on constrain of every ray that will be rendered, so should use SupEpand as mask |
---|
189 | |
---|
190 | mask = Sup ==i; % Not include the Ray that will be use to expand the NonSky |
---|
191 | [yt x] = find(mask); |
---|
192 | |
---|
193 | % find the center point of the sup |
---|
194 | CenterX = round(median(x)); |
---|
195 | CenterY = round(median(yt)); |
---|
196 | |
---|
197 | % if CenterY is lower than ceiling, then set YPointer for that sup to 1 so that it will not come back like a ceiling |
---|
198 | YPointer = [YPointer; CenterY >= ceiling]; |
---|
199 | |
---|
200 | % the horizontal and vertical enforcing in the second step OPT will depends on the height of the sup, so keep record of CenterY |
---|
201 | YPosition = [YPosition; CenterY]; |
---|
202 | |
---|
203 | % Not building PosiM and VarM, have to get rid of NaN depths |
---|
204 | mask(isnan(CleanedDepthMap)) = false; |
---|
205 | SupNuPatch(i) = sum(mask(:)); |
---|
206 | |
---|
207 | % find the center point of the sup |
---|
208 | [yt x] = find(mask); % notice VarM depends on the position of mask if gravity == 1 |
---|
209 | |
---|
210 | if ~all(mask(:)==0) |
---|
211 | if gravity |
---|
212 | if any(CleanedDepthMap(mask) <=0) |
---|
213 | CleanedDepthMap(mask) |
---|
214 | end |
---|
215 | PosiM = blkdiag(PosiM,Posi3D(:,mask)');%.*repmat( mapVert(yt)',[1 3]).*repmat( mapHori(x)',[1 3])); |
---|
216 | VarM = [VarM; VarMap(mask).*(mapVert(yt)').*( mapHori(x)')]; |
---|
217 | else |
---|
218 | PosiM = blkdiag(PosiM,Posi3D(:,mask)'); |
---|
219 | VarM = [VarM; VarMap(mask)]; |
---|
220 | end |
---|
221 | else |
---|
222 | PosiM = blkdiag(PosiM, Posi3D(:,mask)'); |
---|
223 | VarM = [VarM; VarMap(mask)]; |
---|
224 | end |
---|
225 | end |
---|
226 | % set 0 to -1 since this will make the sup higher than the ceiling enforce to come back, if you do not want this just keep 0 as 0 |
---|
227 | YPointer(YPointer==0) = -1; |
---|
228 | |
---|
229 | % buliding CoPlane Matrix========================================================================= |
---|
230 | |
---|
231 | % ===============please ignore here since CoPSTList = [], =============== |
---|
232 | NuSTList = 0; |
---|
233 | if CoPST |
---|
234 | NuSTList = size(CoPSTList,1); |
---|
235 | if ~isempty(CoPSTList) |
---|
236 | [V H] = size(SupNeighborTable); |
---|
237 | SupNeighborTable( CoPSTList(:,1)*V + CoPSTList(:,2)) = 1; |
---|
238 | SupNeighborTable( CoPSTList(:,2)*V + CoPSTList(:,1)) = 1; |
---|
239 | end |
---|
240 | end |
---|
241 | % ================please ignore============================================== |
---|
242 | |
---|
243 | CoPM1 = sparse(0,3*NuSupSize); |
---|
244 | CoPM2 = sparse(0,3*NuSupSize); |
---|
245 | CoPEstDepth = sparse(0,0); |
---|
246 | CoPNorM = []; |
---|
247 | WeiCoP = []; |
---|
248 | if ExtractRelationInfo == 1 |
---|
249 | % keeps the Wei of the relational Coplanar term here |
---|
250 | WeiM = sparse(max(NuSup),max(NuSup)); |
---|
251 | end |
---|
252 | |
---|
253 | for i = NuSup |
---|
254 | % pick the i's neightbot using SupNeighborTable |
---|
255 | Neighbor = find( SupNeighborTable(i,:) ~=0); |
---|
256 | Neighbor = Neighbor( Neighbor> i); |
---|
257 | % setup the relation between sup_i and sup_j since they are neighbors |
---|
258 | for j = Neighbor |
---|
259 | mask = Sup == i; |
---|
260 | SizeMaskAll = sum(mask(:)); |
---|
261 | [y x] = find(mask); |
---|
262 | % Coplanar term only use the center ray of each sup to enfore the cost |
---|
263 | CenterX = round(median(x)); |
---|
264 | CenterY = round(median(y)); |
---|
265 | y = find(mask(:,CenterX)); |
---|
266 | if ~isempty(y) |
---|
267 | CenterY = round(median(y)); |
---|
268 | end |
---|
269 | temp1 = sparse(1, 3*NuSupSize); |
---|
270 | temp2 = sparse(1, 3*NuSupSize); |
---|
271 | temp1(:,(Sup2Para( i)*3-2): Sup2Para( i)*3) = ... |
---|
272 | RayOri(:,CenterY,CenterX)'; |
---|
273 | temp2(:,(Sup2Para( j)*3-2): Sup2Para( j)*3) = ... |
---|
274 | RayOri(:,CenterY,CenterX)'; |
---|
275 | |
---|
276 | % assign wei for each pairs of neighbors |
---|
277 | if MultiScaleFlag |
---|
278 | vector = (MultiScaleSupTable(Sup2Para( i),2:end) == MultiScaleSupTable(Sup2Para( j),2:end)); |
---|
279 | expV = exp(-10*(WeiV*vector' + ShiftCoP) ); |
---|
280 | wei = 1/(1+expV); |
---|
281 | else |
---|
282 | wei = 1; |
---|
283 | end |
---|
284 | |
---|
285 | if ExtractRelationInfo == 1; % keep record |
---|
286 | WeiM(i,j) = wei; |
---|
287 | end |
---|
288 | |
---|
289 | oneRay1 = temp1*wei; |
---|
290 | oneRay2 = temp2*wei; |
---|
291 | |
---|
292 | tempWeiCoP = [SizeMaskAll]; |
---|
293 | CoPEstDepth = [CoPEstDepth; max(median(CleanedDepthMap(mask)),ClosestDist)]; |
---|
294 | |
---|
295 | mask = Sup == j; |
---|
296 | SizeMaskAll = sum(mask(:)); |
---|
297 | [y x] = find(mask); |
---|
298 | |
---|
299 | % Coplanar term only use the center ray of each sup to enfore the cost |
---|
300 | CenterX = round(median(x)); |
---|
301 | CenterY = round(median(y)); |
---|
302 | y = find(mask(:,CenterX)); |
---|
303 | if ~isempty(y) |
---|
304 | CenterY = round(median(y)); |
---|
305 | end |
---|
306 | |
---|
307 | temp1 = sparse(1, 3*NuSupSize); |
---|
308 | temp2 = sparse(1, 3*NuSupSize); |
---|
309 | temp1(:,(Sup2Para( i)*3-2): Sup2Para( i)*3) = ... |
---|
310 | RayOri(:,CenterY,CenterX)'; |
---|
311 | temp2(:,(Sup2Para( j)*3-2): Sup2Para( j)*3) = ... |
---|
312 | RayOri(:,CenterY,CenterX)'; |
---|
313 | |
---|
314 | % Instead of having separate L-1 terms for symmetric co-planar constraint; do the following: |
---|
315 | % If the penaly was ||.||_2^2 + ||.||_2^2; then the co-efficients are |
---|
316 | % some kind of average of two rays. For one norm; we take its average. |
---|
317 | % (do not divide by 2 because the penalty should be double. |
---|
318 | |
---|
319 | CoPM1 = [CoPM1; temp1*wei + oneRay1 ]; |
---|
320 | CoPM2 = [CoPM2; temp2*wei + oneRay2 ]; |
---|
321 | |
---|
322 | tempWeiCoP = [tempWeiCoP; SizeMaskAll]; |
---|
323 | WeiCoP = [WeiCoP; tempWeiCoP]; |
---|
324 | CoPEstDepth = [CoPEstDepth; max(median(CleanedDepthMap(mask)),ClosestDist)]; |
---|
325 | end |
---|
326 | end%========================================================================================================= |
---|
327 | |
---|
328 | %== find the boundary point that might need to be stick ot each other========================================== |
---|
329 | HoriStickM_i = sparse(0,3*NuSupSize); |
---|
330 | HoriStickM_j = sparse(0,3*NuSupSize); |
---|
331 | HoriStickPointInd = []; |
---|
332 | EstDepHoriStick = []; |
---|
333 | |
---|
334 | MAX_POINTS_STITCH_HORI = 2; % the actual code will be modified in another copy of the file |
---|
335 | MIN_POINTS_STITCH = 2; % ERROR: not used. |
---|
336 | |
---|
337 | % ================================================================ |
---|
338 | % NOTE: The actual algorithm should be picking precisely 2 points which are |
---|
339 | % FARTHEST away from the candidate set of neighbors. This algorithm |
---|
340 | % will ALWAYS work and produce no surprising results. |
---|
341 | |
---|
342 | % In some cases, one may experiment with picking only 1 point when the 2 |
---|
343 | % points are too close ---- this will make the algorithm faster; but might |
---|
344 | % produce surprising (e.g. a triangle sticking out) sometimes. |
---|
345 | % An ideal algorithm will reduce the number of points by checking for loops |
---|
346 | % passing through 3 or less superpixels through this matrix; and removing |
---|
347 | % them such that the smallest loop passes through 4 superpixels. (see EE263 |
---|
348 | % for a quick algorithm to do this -- involves product of matrices. |
---|
349 | % ================================================================== |
---|
350 | |
---|
351 | DIST_STICHING_THRESHOLD_HORI = 0.4; |
---|
352 | DIST_STICHING_THRESHOLD_HORI_ONLYCOL = -0.5; % effectively not used, |
---|
353 | |
---|
354 | % Ashutosh added to reduce linear dependancy of the objective ======= |
---|
355 | SupPixelNeighborList = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
356 | SupPixelParsedList = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
357 | recordAdded1 = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
358 | recordAdded2 = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
359 | addedIndexList = [ ]; |
---|
360 | % =================================================================== |
---|
361 | |
---|
362 | % locate sup boundary at horizontal and vertical dircetion only for stitching terms |
---|
363 | BounaryPHori = conv2(Sup,[1 -1],'same') ~=0; |
---|
364 | BounaryPHori(:,end) = 0; |
---|
365 | BounaryPVert = conv2(Sup,[1; -1],'same') ~=0; |
---|
366 | BounaryPVert(end,:) = 0; |
---|
367 | |
---|
368 | % boundariesHoriIndex = find(BounaryPHori==1)'; |
---|
369 | % pre-select the boundary in order with the NuSup order |
---|
370 | for l = NuSup |
---|
371 | mask = Sup == l; |
---|
372 | boundariesHoriIndex = find(BounaryPHori==1 & mask)'; |
---|
373 | for i = boundariesHoriIndex |
---|
374 | j = i+Default.VertYNuDepth; |
---|
375 | if Sup(i) == 0 || Sup(j) == 0 |
---|
376 | continue; |
---|
377 | end |
---|
378 | SupPixelParsedList(Sup(i),Sup(j)) = SupPixelParsedList(Sup(i),Sup(j)) + 1; |
---|
379 | |
---|
380 | if SupPixelNeighborList(Sup(i),Sup(j)) == 0 |
---|
381 | recordAdded1(Sup(i),Sup(j)) = i; |
---|
382 | elseif SupPixelNeighborList(Sup(i),Sup(j)) >= MAX_POINTS_STITCH_HORI |
---|
383 | continue; |
---|
384 | elseif SupPixelNeighborList(Sup(i),Sup(j)) == 1 % inside this remove the close stiching terms |
---|
385 | rowN = ceil(i/55); colN = rem(i,55); |
---|
386 | rowN_older = ceil( recordAdded1(Sup(i),Sup(j)) / 55); |
---|
387 | colN_older = rem( recordAdded1(Sup(i),Sup(j)), 55); |
---|
388 | if abs(rowN - rowN_older) + (55/305)*abs(colN - colN_older) > DIST_STICHING_THRESHOLD_HORI && ... |
---|
389 | abs(colN - colN_older) > DIST_STICHING_THRESHOLD_HORI_ONLYCOL |
---|
390 | recordAdded2(Sup(i),Sup(j)) = i; |
---|
391 | else |
---|
392 | continue; |
---|
393 | end |
---|
394 | elseif SupPixelNeighborList(Sup(i),Sup(j)) == 2 %Assuming MAX_POINTS_STITCH = 3 |
---|
395 | rowN = ceil(i/55); colN = rem(i,55); |
---|
396 | rowN_older1 = ceil( recordAdded1(Sup(i),Sup(j)) / 55); |
---|
397 | colN_older1 = rem( recordAdded1(Sup(i),Sup(j)), 55); |
---|
398 | rowN_older2 = ceil( recordAdded2(Sup(i),Sup(j)) / 55); |
---|
399 | colN_older2 = rem( recordAdded2(Sup(i),Sup(j)), 55); |
---|
400 | |
---|
401 | if abs(rowN - rowN_older1) + (55/305)*abs(colN - colN_older1) > DIST_STICHING_THRESHOLD_HORI && ... |
---|
402 | abs(rowN - rowN_older2) + (55/305)*abs(colN - colN_older2) > DIST_STICHING_THRESHOLD_HORI |
---|
403 | ; |
---|
404 | else |
---|
405 | continue; |
---|
406 | end |
---|
407 | end |
---|
408 | |
---|
409 | % If you come here, it means you are probably adding it. |
---|
410 | SupPixelNeighborList(Sup(i),Sup(j)) = SupPixelNeighborList(Sup(i),Sup(j)) + 1; |
---|
411 | addedIndexList = [addedIndexList i]; |
---|
412 | end |
---|
413 | end |
---|
414 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
415 | %%%% this is supposed to be where the y and nu information are computed, |
---|
416 | %%%% keep an eye on this. |
---|
417 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
418 | % after analysis finally add to HoriStickM_j |
---|
419 | WeightHoriNeighborStitch = [ ]; |
---|
420 | for i = addedIndexList |
---|
421 | j = i+Default.VertYNuDepth; |
---|
422 | WeightHoriNeighborStitch = [WeightHoriNeighborStitch; SupPixelParsedList(Sup(i),Sup(j)) / ... |
---|
423 | SupPixelNeighborList(Sup(i),Sup(j)) ]; |
---|
424 | |
---|
425 | Target(1) = Sup2Para(Sup(i)); |
---|
426 | Target(2) = Sup2Para(Sup(j)); |
---|
427 | rayBoundary(:,1) = RayOri(:,i); |
---|
428 | rayBoundary(:,2) = RayOri(:,i); |
---|
429 | if MultiScaleFlag |
---|
430 | vector = (MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end)); |
---|
431 | expV = exp(-10*(WeiV*vector' + ShiftStick) ); |
---|
432 | betaTemp = StickHori*(0.5+1/(1+expV)); |
---|
433 | % therr should always be sticking (know don't care about occlusion) |
---|
434 | else |
---|
435 | betaTemp = StickHori; |
---|
436 | end |
---|
437 | temp = sparse(3,NuSupSize); |
---|
438 | temp(:,Target(1)) = rayBoundary(:,1); |
---|
439 | HoriStickM_i = [HoriStickM_i; betaTemp*temp(:)']; |
---|
440 | temp = sparse(3,NuSupSize); |
---|
441 | temp(:,Target(2)) = rayBoundary(:,2); |
---|
442 | HoriStickM_j = [HoriStickM_j; betaTemp*temp(:)']; |
---|
443 | EstDepHoriStick = [EstDepHoriStick; sqrt(max(CleanedDepthMap(i),ClosestDist)*max(CleanedDepthMap(j),ClosestDist))]; |
---|
444 | HoriStickPointInd = [HoriStickPointInd i ]; |
---|
445 | end |
---|
446 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
447 | %%%% end of y and nu computation. |
---|
448 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
449 | % ============================================== |
---|
450 | |
---|
451 | % ======== finding the unucessary stiching points in Vertical direction ==== |
---|
452 | VertStickM_i = sparse(0,3*NuSupSize); |
---|
453 | VertStickM_j = sparse(0,3*NuSupSize); |
---|
454 | VertStickPointInd = []; |
---|
455 | EstDepVertStick = []; |
---|
456 | |
---|
457 | MAX_POINTS_STITCH_VERT = 4; %3 |
---|
458 | DIST_STICHING_THRESHOLD_VERT = 0.1; %0.3 |
---|
459 | DIST_STICHING_THRESHOLD_VERT_ONLYCOL = -0.5; % effectively not used, ideally should be 0.5; i.e., the point should be farther in col direction because that is the direction of the edge. |
---|
460 | |
---|
461 | % Ashutosh added to reduce linear dependancy of the objective ======= |
---|
462 | SupPixelNeighborList = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
463 | SupPixelParsedList = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
464 | recordAdded1 = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
465 | recordAdded2 = sparse( max(Sup(:)), max(Sup(:)) ); |
---|
466 | addedIndexList = [ ]; |
---|
467 | % =============================== |
---|
468 | |
---|
469 | % pre-select the boundary in order with the NuSup order |
---|
470 | for l = NuSup |
---|
471 | mask = Sup == l; |
---|
472 | for i = find(BounaryPVert==1 & mask)' |
---|
473 | j = i+1; |
---|
474 | if Sup(i) == 0 || Sup(j) == 0 |
---|
475 | continue; |
---|
476 | end |
---|
477 | SupPixelParsedList(Sup(i),Sup(j)) = SupPixelParsedList(Sup(i),Sup(j)) + 1; |
---|
478 | |
---|
479 | if SupPixelNeighborList(Sup(i),Sup(j)) == 0 |
---|
480 | recordAdded1(Sup(i),Sup(j)) = i; |
---|
481 | elseif SupPixelNeighborList(Sup(i),Sup(j)) >= MAX_POINTS_STITCH_VERT |
---|
482 | continue; |
---|
483 | elseif SupPixelNeighborList(Sup(i),Sup(j)) == 1 % inside this remove the close stiching terms |
---|
484 | rowN = ceil(i/55); colN = rem(i,55); |
---|
485 | rowN_older = ceil( recordAdded1(Sup(i),Sup(j)) / 55); |
---|
486 | colN_older = rem( recordAdded1(Sup(i),Sup(j)), 55); |
---|
487 | if abs(rowN - rowN_older) + (55/305)*abs(colN - colN_older) > DIST_STICHING_THRESHOLD_VERT && ... |
---|
488 | abs(colN - colN_older) > DIST_STICHING_THRESHOLD_VERT_ONLYCOL |
---|
489 | recordAdded2(Sup(i),Sup(j)) = i; |
---|
490 | else |
---|
491 | continue; |
---|
492 | end |
---|
493 | elseif SupPixelNeighborList(Sup(i),Sup(j)) == 2 %Assuming MAX_POINTS_STITCH = 3 |
---|
494 | rowN = ceil(i/55); colN = rem(i,55); |
---|
495 | rowN_older1 = ceil( recordAdded1(Sup(i),Sup(j)) / 55); |
---|
496 | colN_older1 = rem( recordAdded1(Sup(i),Sup(j)), 55); |
---|
497 | rowN_older2 = ceil( recordAdded2(Sup(i),Sup(j)) / 55); |
---|
498 | colN_older2 = rem( recordAdded2(Sup(i),Sup(j)), 55); |
---|
499 | |
---|
500 | if abs(rowN - rowN_older1) + (55/305)*abs(colN - colN_older1) > DIST_STICHING_THRESHOLD_VERT && ... |
---|
501 | abs(rowN - rowN_older2) + (55/305)*abs(colN - colN_older2) > DIST_STICHING_THRESHOLD_VERT |
---|
502 | ; |
---|
503 | else |
---|
504 | continue; |
---|
505 | end |
---|
506 | end |
---|
507 | |
---|
508 | % If you come here, it means you are probably adding it. |
---|
509 | SupPixelNeighborList(Sup(i),Sup(j)) = SupPixelNeighborList(Sup(i),Sup(j)) + 1; |
---|
510 | addedIndexList = [addedIndexList i]; |
---|
511 | end |
---|
512 | end |
---|
513 | |
---|
514 | % after analysis finally add to VertStickM_j |
---|
515 | WeightVertNeighborStitch = [ ]; |
---|
516 | for i = addedIndexList |
---|
517 | j = i+1; |
---|
518 | WeightVertNeighborStitch = [WeightVertNeighborStitch; SupPixelParsedList(Sup(i),Sup(j)) / ... |
---|
519 | SupPixelNeighborList(Sup(i),Sup(j)) ]; |
---|
520 | Target(1) = Sup2Para(Sup(i)); |
---|
521 | Target(2) = Sup2Para(Sup(j)); |
---|
522 | rayBoundary(:,1) = RayOri(:,i); |
---|
523 | rayBoundary(:,2) = RayOri(:,i); |
---|
524 | if MultiScaleFlag |
---|
525 | vector = (MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end)); |
---|
526 | expV = exp(-10*(WeiV*vector' + ShiftStick) ); |
---|
527 | betaTemp = StickVert*(0.5+1/(1+expV)); |
---|
528 | % therr should always be sticking (know don't care about occlusion) |
---|
529 | else |
---|
530 | betaTemp = StickVert; |
---|
531 | end |
---|
532 | temp = sparse(3,NuSupSize); |
---|
533 | temp(:,Target(1)) = rayBoundary(:,1); |
---|
534 | VertStickM_i = [VertStickM_i; betaTemp*temp(:)']; |
---|
535 | temp = sparse(3,NuSupSize); |
---|
536 | temp(:,Target(2)) = rayBoundary(:,2); |
---|
537 | VertStickM_j = [VertStickM_j; betaTemp*temp(:)']; |
---|
538 | EstDepVertStick = [EstDepVertStick; sqrt(max(CleanedDepthMap(i),ClosestDist)*max(CleanedDepthMap(j),ClosestDist))]; |
---|
539 | VertStickPointInd = [VertStickPointInd i ]; |
---|
540 | end |
---|
541 | % ====finished finding the unucessary stiching points in Vertical direction === |
---|
542 | |
---|
543 | % ======================================Finish building up matrix=====================hard work====================== |
---|
544 | |
---|
545 | |
---|
546 | % ================================================================================================================ |
---|
547 | % Start Decompose the image align with superpixels ================================================================================== |
---|
548 | % define the decomposition in both X and Y direction |
---|
549 | |
---|
550 | TotalRectX = 2*XNuDecompose-1; |
---|
551 | TotalRectY= 2*YNuDecompose-1; |
---|
552 | PlanePara = NaN*ones(3*NuSupSize,1); % setup the lookup table for the solved plane parameter |
---|
553 | |
---|
554 | for k = 0:(TotalRectX-1) |
---|
555 | l = rem(k*2,(TotalRectX)); |
---|
556 | RangeX = (1+ceil(Default.HoriXNuDepth/XNuDecompose)*l/2):... |
---|
557 | min((1+ceil(Default.HoriXNuDepth/XNuDecompose)*(l/2+1)),Default.HoriXNuDepth); |
---|
558 | RangeX = ceil(RangeX); |
---|
559 | for q = 0:(TotalRectY-1) |
---|
560 | l = rem(q*2,(TotalRectY)); |
---|
561 | RangeY = (1+ceil(Default.VertYNuDepth/YNuDecompose)*l/2):... |
---|
562 | min((1+ceil(Default.VertYNuDepth/YNuDecompose)*(l/2+1)),Default.VertYNuDepth); |
---|
563 | RangeY = ceil(RangeY); |
---|
564 | mask = zeros(size(Sup)); |
---|
565 | mask(RangeY,RangeX) = 1; |
---|
566 | mask =logical(mask); |
---|
567 | SubSup = sort(setdiff(unique( reshape( Sup(RangeY,RangeX),1,[])),0)); |
---|
568 | BoundarySup = []; |
---|
569 | |
---|
570 | if any(SubSup <=0) |
---|
571 | SubSup(SubSup<=0) |
---|
572 | end |
---|
573 | BoundarySup = find(sum(SupNeighborTable(SubSup,:), 1) ~=0); |
---|
574 | BoundarySup = unique(setdiff(BoundarySup,[0 SubSup] )); |
---|
575 | |
---|
576 | % chech if BoundarySup non-NaN in PlanePara |
---|
577 | checkNoNNaN = ~isnan(PlanePara(Sup2Para(BoundarySup)*3)); |
---|
578 | BoundarySup = BoundarySup(checkNoNNaN); |
---|
579 | TotalSup = sort([SubSup BoundarySup]); |
---|
580 | |
---|
581 | % define the pointer to extract the data for this specific decomposition |
---|
582 | SubSupPtr = [ Sup2Para(SubSup)*3-2;... |
---|
583 | Sup2Para(SubSup)*3-1;... |
---|
584 | Sup2Para(SubSup)*3]; |
---|
585 | SubSupPtr = SubSupPtr(:); |
---|
586 | BoundarySupPtr = [ Sup2Para(BoundarySup)*3-2;... |
---|
587 | Sup2Para(BoundarySup)*3-1;... |
---|
588 | Sup2Para(BoundarySup)*3]; |
---|
589 | BoundarySupPtr = BoundarySupPtr(:); |
---|
590 | NuSubSupSize = size(SubSup,2); |
---|
591 | |
---|
592 | % simply extract NewRayAllM NewPosiM NewCoPM ========================================= |
---|
593 | % NewHoriStickM NewVertStickM for the specific decomposition |
---|
594 | NewRayAllM = RayAllM(:,SubSupPtr); |
---|
595 | tar = sum(NewRayAllM ~= 0,2) == 3; |
---|
596 | NewRayAllM = NewRayAllM(tar,:); |
---|
597 | NewPosiM = PosiM(:,SubSupPtr); |
---|
598 | tar = sum(NewPosiM ~= 0,2) == 3; |
---|
599 | NewPosiM = NewPosiM(tar,:); |
---|
600 | NewVarM = VarM(tar); |
---|
601 | |
---|
602 | NewCoPM = CoPM1(:,SubSupPtr) - CoPM2(:,SubSupPtr); |
---|
603 | NewCoPMBound = CoPM1(:,BoundarySupPtr) - CoPM2(:,BoundarySupPtr); |
---|
604 | |
---|
605 | tar = sum(NewCoPM ~= 0,2) + sum(NewCoPMBound ~= 0,2)==6; |
---|
606 | NewCoPMBound = NewCoPMBound*PlanePara(BoundarySupPtr); |
---|
607 | NewCoPM = NewCoPM(tar,:); |
---|
608 | NewCoPMBound = NewCoPMBound(tar); |
---|
609 | NewCoPEstDepth = CoPEstDepth(tar); |
---|
610 | |
---|
611 | NewHoriStickM = HoriStickM_i(:,SubSupPtr)-HoriStickM_j(:,SubSupPtr); |
---|
612 | NewHoriStickMBound = HoriStickM_i(:,BoundarySupPtr)-HoriStickM_j(:,BoundarySupPtr); |
---|
613 | |
---|
614 | tar = sum(NewHoriStickM ~= 0,2)+ sum(NewHoriStickMBound ~= 0,2) ==6; |
---|
615 | NewHoriStickM = NewHoriStickM(tar,:); |
---|
616 | NewEstDepHoriStick = EstDepHoriStick(tar); |
---|
617 | NewHoriStickMBound = NewHoriStickMBound*PlanePara(BoundarySupPtr); |
---|
618 | NewHoriStickMBound = NewHoriStickMBound(tar); |
---|
619 | NewWeightHoriNeighborStitch = WeightHoriNeighborStitch(tar); |
---|
620 | |
---|
621 | NewVertStickM = VertStickM_i(:,SubSupPtr)-VertStickM_j(:,SubSupPtr); |
---|
622 | NewVertStickMBound = VertStickM_i(:,BoundarySupPtr)-VertStickM_j(:,BoundarySupPtr); |
---|
623 | |
---|
624 | tar = sum(NewVertStickM ~= 0,2) + sum(NewVertStickMBound ~= 0,2)==6; |
---|
625 | NewVertStickM = NewVertStickM(tar,:); |
---|
626 | NewEstDepVertStick = EstDepVertStick(tar); |
---|
627 | NewVertStickMBound = NewVertStickMBound*PlanePara(BoundarySupPtr); |
---|
628 | NewVertStickMBound = NewVertStickMBound(tar); |
---|
629 | NewWeightVertNeighborStitch = WeightVertNeighborStitch(tar); |
---|
630 | |
---|
631 | WeightsSelfTerm = 1 ./ exp(abs(NewVarM)/BandWith); |
---|
632 | % ==========================end of extraction ==================== |
---|
633 | |
---|
634 | fprintf([' ' num2str( toc(inferenceTime) ) '\n : In 1st level Optimization, using new solver.' ... |
---|
635 | '(' num2str(k+1) '/' num2str((TotalRectX-1)+1) ',' num2str(l+1) '/' num2str((TotalRectY-1)+1) ')']); |
---|
636 | |
---|
637 | global A b S inq |
---|
638 | |
---|
639 | A = [ sparse(1:length(WeightsSelfTerm),1:length(WeightsSelfTerm),WeightsSelfTerm ) * NewPosiM;... |
---|
640 | sparse(1:length(NewCoPEstDepth),1:length(NewCoPEstDepth),NewCoPEstDepth*Center ) * NewCoPM;... |
---|
641 | sparse(1:length(NewEstDepHoriStick), 1:length(NewEstDepHoriStick), NewEstDepHoriStick.*NewWeightHoriNeighborStitch) * NewHoriStickM;... |
---|
642 | sparse(1:length(NewEstDepVertStick), 1:length(NewEstDepVertStick), NewEstDepVertStick.*NewWeightVertNeighborStitch) * NewVertStickM;... |
---|
643 | ]; |
---|
644 | |
---|
645 | b = [ ones(size(NewPosiM,1),1) .* WeightsSelfTerm;... |
---|
646 | -NewCoPMBound.*NewCoPEstDepth*Center;... |
---|
647 | -NewHoriStickMBound.*NewEstDepHoriStick.*NewWeightHoriNeighborStitch;... |
---|
648 | -NewVertStickMBound.*NewEstDepVertStick.*NewWeightVertNeighborStitch]; |
---|
649 | |
---|
650 | temp = zeros(1, NuSubSupSize*3); |
---|
651 | temp(3*(1:NuSubSupSize)-1) = YPointer(Sup2Para(SubSup)); |
---|
652 | temp = sparse(1:length(temp), 1:length(temp), temp); |
---|
653 | temp( sum(temp,2) ==0,:) = []; |
---|
654 | S = [temp;... |
---|
655 | NewRayAllM;... |
---|
656 | -NewRayAllM]; |
---|
657 | inq = [ sparse(size(temp,1), 1);... |
---|
658 | - 1/ClosestDist*ones(size(NewRayAllM,1),1);... |
---|
659 | 1/FarestDist*ones(size(NewRayAllM,1),1)]; |
---|
660 | Para.ClosestDist = ClosestDist; |
---|
661 | Para.FarestDist = FarestDist; |
---|
662 | Para.ptry = []; |
---|
663 | Para.ptrz = []; |
---|
664 | Para.Dist_Start = []; |
---|
665 | |
---|
666 | % solve the OPT problem using the new solver |
---|
667 | [x_ashIterator, alfa, status, history, T_nt_hist] = ... |
---|
668 | SigmoidLogBarrierSolver( Para, [], [], [], '', [], [], solverVerboseLevel); |
---|
669 | |
---|
670 | % check if the constrian still satisfied |
---|
671 | if any(S*x_ashIterator+inq > 0 ) |
---|
672 | disp('Inequality not satisfied'); |
---|
673 | max( S*x_ashIterator+inq) |
---|
674 | elseif status == 2 |
---|
675 | fprintf([' Success with alfa=' num2str(alfa)]); |
---|
676 | end |
---|
677 | |
---|
678 | % assign the solution of specific decomposition to the corresponding entries of the whole problems |
---|
679 | PlanePara(SubSupPtr) = x_ashIterator; |
---|
680 | end |
---|
681 | end |
---|
682 | |
---|
683 | % build the whole image |
---|
684 | PlanePara = reshape(PlanePara,3,[]); |
---|
685 | % porject the ray on planes to generate the ProjDepth |
---|
686 | FitDepthPPCP = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth); |
---|
687 | FitDepthPPCP(~maskSkyEroded) = (1./sum(PlanePara(:,Sup2Para(SupEpand(~maskSkyEroded ))).*Ray(:,~maskSkyEroded ),1))'; |
---|
688 | FitDepthPPCP = reshape(FitDepthPPCP,Default.VertYNuDepth,[]); |
---|
689 | [Position3DFitedPPCP] = im_cr2w_cr(FitDepthPPCP,permute(Ray,[2 3 1])); |
---|
690 | |
---|
691 | if NoSecondStep % if no second step |
---|
692 | Position3DFitedPPCP(3,:) = -Position3DFitedPPCP(3,:); |
---|
693 | Position3DFitedPPCP = permute(Position3DFitedPPCP,[2 3 1]); |
---|
694 | RR =permute(Ray,[2 3 1]); |
---|
695 | temp = RR(:,:,1:2)./repmat(RR(:,:,3),[1 1 2]); |
---|
696 | PositionTex = permute(temp./repmat(cat(3,Default.a_default,Default.b_default),[Default.VertYNuDepth Default.HoriXNuDepth 1])+repmat(cat(3,Default.Ox_default,Default.Oy_default),[Default.VertYNuDepth Default.HoriXNuDepth 1]),[3 1 2]); |
---|
697 | PositionTex = permute(PositionTex,[2 3 1]) |
---|
698 | disp('First level Wrl'); |
---|
699 | |
---|
700 | % write wrl file |
---|
701 | WrlFacestHroiReduce(Position3DFitedPPCP,PositionTex,SupOri, [ Default.filename{1} '1st'],[ Default.filename{1} '1st'], ... |
---|
702 | Default.OutPutFolder, GridFlag, 0); |
---|
703 | system(['gzip -9 -c ' Default.OutPutFolder Default.filename{1} '1st.wrl > ' ... |
---|
704 | Default.OutPutFolder Default.filename{1} '1st.wrl.gz']); |
---|
705 | copyfile([ Default.OutPutFolder Default.filename{1} '1st.wrl.gz '], ... |
---|
706 | [ Default.OutPutFolder Default.filename{1} '1st.wrl'],'f'); |
---|
707 | delete([Default.OutPutFolder Default.filename{1} '1st.wrl.gz']); |
---|
708 | end |
---|
709 | %==================Finished for one step MRF================================================================================================== |
---|
710 | |
---|
711 | if NoSecondStep |
---|
712 | return; |
---|
713 | end |
---|
714 | |
---|
715 | % ====================================following are 2nd step MRF to give more visually pleasing result======================================= |
---|
716 | % generating new PosiMPPCP using the new position |
---|
717 | normPara = norms(PlanePara); |
---|
718 | normalizedPara = PlanePara ./ repmat( normPara, [3 1]); |
---|
719 | groundPara = abs(normalizedPara(2,:)) >= groundThreshold(YPosition); |
---|
720 | groundParaInd = find(groundPara); |
---|
721 | verticalPara = abs(normalizedPara(2,:)) <= verticalThreshold(YPosition); % change to have different range of vertical thre ================ |
---|
722 | verticalParaInd = find(verticalPara); |
---|
723 | |
---|
724 | indexVertical = find( verticalPara)*3-1; |
---|
725 | indexGroundX = find( groundPara)*3-2; |
---|
726 | indexGroundZ = find( groundPara)*3; |
---|
727 | |
---|
728 | PosiMPPCP = sparse(0,0); |
---|
729 | VarM2 = sparse(0,0); |
---|
730 | |
---|
731 | % forming new supporting matrix using new depth and get rid of the support of the vertical plane |
---|
732 | for i = NuSup |
---|
733 | mask = Sup == i; |
---|
734 | if any(verticalParaInd == Sup2Para(i)) & FlagRemoveVerticalSupport % Remove Vertical Depths Supports |
---|
735 | mask = logical(zeros(size(Sup))); |
---|
736 | end |
---|
737 | PosiMPPCP = blkdiag(PosiMPPCP, Position3DFitedPPCP(:,mask)'); |
---|
738 | VarM2 = [VarM2; VarMap(mask)]; |
---|
739 | if length( find(isnan(PosiMPPCP)) ) |
---|
740 | disp('PosiMPPCP is NaN'); |
---|
741 | end |
---|
742 | end |
---|
743 | |
---|
744 | % Start Decompose image ========================= |
---|
745 | TotalRectX = 2*XNuDecompose-1; |
---|
746 | TotalRectY = 2*YNuDecompose-1; |
---|
747 | PlanePara = NaN*ones(3*NuSupSize,1); % setup the lookuptable for the solved plane parameter |
---|
748 | |
---|
749 | for k = 0:(TotalRectX-1) |
---|
750 | l = rem(k*2,(TotalRectX)); |
---|
751 | RangeX = (1+ceil(Default.HoriXNuDepth/XNuDecompose)*l/2):... |
---|
752 | min((1+ceil(Default.HoriXNuDepth/XNuDecompose)*(l/2+1)),Default.HoriXNuDepth); |
---|
753 | RangeX = ceil(RangeX); |
---|
754 | for q = 0:(TotalRectY-1) |
---|
755 | l = rem(q*2,(TotalRectY)); |
---|
756 | RangeY = (1+ceil(Default.VertYNuDepth/YNuDecompose)*l/2):... |
---|
757 | min((1+ceil(Default.VertYNuDepth/YNuDecompose)*(l/2+1)),Default.VertYNuDepth); |
---|
758 | RangeY = ceil(RangeY); |
---|
759 | mask = zeros(size(Sup)); |
---|
760 | mask(RangeY,RangeX) = 1; |
---|
761 | mask =logical(mask); |
---|
762 | SubSup = sort(setdiff( unique( reshape( Sup(RangeY,RangeX),1,[])),0)); |
---|
763 | BoundarySup = []; |
---|
764 | BoundarySup = find(sum(SupNeighborTable(SubSup,:), 1) ~=0); |
---|
765 | BoundarySup = unique(setdiff(BoundarySup,[0 SubSup] )); |
---|
766 | |
---|
767 | % chech if BoundarySup non-NaN in PlanePara |
---|
768 | checkNoNNaN = ~isnan(PlanePara(Sup2Para(BoundarySup)*3)); |
---|
769 | BoundarySup = BoundarySup(checkNoNNaN); |
---|
770 | TotalSup = sort([SubSup BoundarySup]); |
---|
771 | |
---|
772 | SubSupPtr = [ Sup2Para(SubSup)*3-2;... |
---|
773 | Sup2Para(SubSup)*3-1;... |
---|
774 | Sup2Para(SubSup)*3]; |
---|
775 | SubSupPtr = SubSupPtr(:); |
---|
776 | BoundarySupPtr = [ Sup2Para(BoundarySup)*3-2;... |
---|
777 | Sup2Para(BoundarySup)*3-1;... |
---|
778 | Sup2Para(BoundarySup)*3]; |
---|
779 | BoundarySupPtr =BoundarySupPtr(:); |
---|
780 | NuSubSupSize = size(SubSup,2); |
---|
781 | SubSup2Para = sparse(1,max(SubSup)); |
---|
782 | SubSup2Para(SubSup) = 1:NuSubSupSize; |
---|
783 | |
---|
784 | % clearn RayAllM PosiM CoPM1 HoriStickM_i VertStickM_i |
---|
785 | NewRayAllM = RayAllM(:,SubSupPtr); |
---|
786 | tar = sum(NewRayAllM ~= 0,2) ==3; |
---|
787 | NewRayAllM = NewRayAllM(tar,:); |
---|
788 | |
---|
789 | NewPosiMPPCP = PosiMPPCP(:,SubSupPtr); |
---|
790 | tar = sum(NewPosiMPPCP ~= 0,2) ==3; |
---|
791 | NewPosiMPPCP = NewPosiMPPCP(tar,:); |
---|
792 | NewVarM = VarM2(tar); |
---|
793 | |
---|
794 | NewCoPM = CoPM1(:,SubSupPtr) - CoPM2(:,SubSupPtr); |
---|
795 | NewCoPMBound = CoPM1(:,BoundarySupPtr) - CoPM2(:,BoundarySupPtr); |
---|
796 | tar = sum( NewCoPM ~= 0,2) + sum( NewCoPMBound ~= 0,2) ==6; |
---|
797 | NewCoPM = NewCoPM(tar,:); |
---|
798 | NewCoPMBound = NewCoPMBound*PlanePara(BoundarySupPtr); % column vertor |
---|
799 | NewCoPMBound = NewCoPMBound(tar); |
---|
800 | NewCoPEstDepth = CoPEstDepth(tar); |
---|
801 | |
---|
802 | NewHoriStickM = HoriStickM_i(:,SubSupPtr)-HoriStickM_j(:,SubSupPtr); |
---|
803 | NewHoriStickMBound = HoriStickM_i(:,BoundarySupPtr)-HoriStickM_j(:,BoundarySupPtr); |
---|
804 | tar = sum(NewHoriStickM ~= 0,2) + sum( NewHoriStickMBound ~= 0,2)==6; |
---|
805 | NewHoriStickM = NewHoriStickM(tar,:); |
---|
806 | NewHoriStickMBound = NewHoriStickMBound*PlanePara(BoundarySupPtr); % column vertor |
---|
807 | NewHoriStickMBound = NewHoriStickMBound(tar); |
---|
808 | NewEstDepHoriStick = EstDepHoriStick(tar); |
---|
809 | NewWeightHoriNeighborStitch = WeightHoriNeighborStitch(tar); |
---|
810 | |
---|
811 | NewVertStickM = VertStickM_i(:,SubSupPtr)-VertStickM_j(:,SubSupPtr); |
---|
812 | NewVertStickMBound = VertStickM_i(:, BoundarySupPtr)-VertStickM_j(:,BoundarySupPtr); |
---|
813 | tar = sum(NewVertStickM ~= 0,2) + sum(NewVertStickMBound ~= 0,2)==6; |
---|
814 | NewVertStickM = NewVertStickM(tar,:); |
---|
815 | NewVertStickMBound = NewVertStickMBound*PlanePara(BoundarySupPtr); % column vertor |
---|
816 | NewVertStickMBound = NewVertStickMBound(tar); |
---|
817 | NewEstDepVertStick = EstDepVertStick(tar); |
---|
818 | NewWeightVertNeighborStitch = WeightVertNeighborStitch(tar); |
---|
819 | |
---|
820 | % try reduce the vertical constrain |
---|
821 | NonVertPtr = setdiff( 1:(3*NuSubSupSize), SubSup2Para( NuSup( (intersect( indexVertical,SubSupPtr ) +1)/3))*3-1 ); |
---|
822 | YNoComingBack = YPointer(Sup2Para(SubSup)); |
---|
823 | YNoComingBack(SubSup2Para( NuSup( (intersect( indexVertical,SubSupPtr ) +1)/3))) = []; |
---|
824 | YCompMask = zeros(1,3*NuSubSupSize); |
---|
825 | YCompMask(3*(1:NuSubSupSize)-1) = 1; |
---|
826 | YCompMask = YCompMask(NonVertPtr); |
---|
827 | XCompMask = zeros(1,3*NuSubSupSize); |
---|
828 | XCompMask( SubSup2Para( NuSup( (intersect( indexGroundX,SubSupPtr ) +2)/3))*3-2 ) = 1; |
---|
829 | XCompMask = XCompMask(NonVertPtr); |
---|
830 | ZCompMask = zeros(1,3*NuSubSupSize); |
---|
831 | ZCompMask( SubSup2Para( NuSup( (intersect( indexGroundZ,SubSupPtr ) )/3))*3 ) = 1; |
---|
832 | ZCompMask = ZCompMask(NonVertPtr); |
---|
833 | GroundMask = intersect( find(groundPara), ( SubSupPtr( 3:3:( size(SubSupPtr,1)) )/3 ) ) ; |
---|
834 | |
---|
835 | % version written by Ashutosh |
---|
836 | |
---|
837 | tempGroundX = sparse(1:length(XCompMask), 1:length(XCompMask), XCompMask); |
---|
838 | tempGroundX( logical(XCompMask) ) = tempGroundX( logical(XCompMask) )./normPara( GroundMask ); |
---|
839 | tempGroundX( sum(tempGroundX,2) == 0,:) = []; |
---|
840 | tempGroundZ = sparse(1:length(ZCompMask), 1:length(ZCompMask), ZCompMask ); |
---|
841 | tempGroundZ( logical(ZCompMask) ) = tempGroundZ( logical(ZCompMask) )./normPara( GroundMask ); |
---|
842 | tempGroundZ( sum(tempGroundZ,2) == 0,:)= []; |
---|
843 | |
---|
844 | A = [ sparse(1:length(NewVarM),1:length(NewVarM),1./exp(abs(NewVarM)/BandWith)) * NewPosiMPPCP(:, NonVertPtr);... |
---|
845 | sparse(1:length(NewCoPEstDepth), 1:length(NewCoPEstDepth), NewCoPEstDepth * Center) * NewCoPM(:, NonVertPtr);... |
---|
846 | sparse(1:length(NewEstDepHoriStick), 1:length(NewEstDepHoriStick), NewEstDepHoriStick.*NewWeightHoriNeighborStitch) * NewHoriStickM(:, NonVertPtr);... |
---|
847 | sparse(1:length(NewEstDepVertStick), 1:length(NewEstDepVertStick), NewEstDepVertStick.*NewWeightVertNeighborStitch) * NewVertStickM(:, NonVertPtr);... |
---|
848 | tempGroundX;... |
---|
849 | tempGroundZ... |
---|
850 | ]; |
---|
851 | % +10*norm( ( Para( logical(XCompMask) ))./... |
---|
852 | % normPara( GroundMask )', 1)... |
---|
853 | % +10*norm( ( Para( logical(ZCompMask)))./... |
---|
854 | % normPara( GroundMask )', 1) ... |
---|
855 | |
---|
856 | b = [ 1./exp(abs(NewVarM)/BandWith); ... |
---|
857 | -Center*NewCoPMBound.*NewCoPEstDepth; ... |
---|
858 | -NewHoriStickMBound.*NewEstDepHoriStick.*NewWeightHoriNeighborStitch;... |
---|
859 | -NewVertStickMBound.*NewEstDepVertStick.*NewWeightVertNeighborStitch;... |
---|
860 | sparse(size(tempGroundX,1),1);... |
---|
861 | sparse(size(tempGroundZ,1),1)... |
---|
862 | ]; |
---|
863 | |
---|
864 | temp = YCompMask; |
---|
865 | temp(logical(YCompMask)) = YNoComingBack; |
---|
866 | temp = sparse(1:length(temp), 1:length(temp), temp); |
---|
867 | temp( sum(temp,2) ==0,:) = []; |
---|
868 | S = [ temp;... |
---|
869 | NewRayAllM(:,NonVertPtr);... |
---|
870 | -NewRayAllM(:,NonVertPtr);... |
---|
871 | ]; |
---|
872 | inq = [ sparse(size(temp,1), 1);... |
---|
873 | - 1/ClosestDist*ones(size(NewRayAllM,1),1);... |
---|
874 | 1/FarestDist*ones(size(NewRayAllM,1),1);... |
---|
875 | ]; |
---|
876 | |
---|
877 | Para.ClosestDist = ClosestDist; |
---|
878 | Para.FarestDist = FarestDist; |
---|
879 | |
---|
880 | % build up ptry and ptrz adapt fot NonVertPtr |
---|
881 | Para.ptry = zeros(size(NewRayAllM,2),1); |
---|
882 | Para.ptry(2:3:size(NewRayAllM,2)) = 1; |
---|
883 | Para.ptry = logical(Para.ptry(NonVertPtr)); |
---|
884 | Para.ptrz = zeros(size(NewRayAllM,2),1); |
---|
885 | Para.ptrz(3:3:size(NewRayAllM,2)) = 1; |
---|
886 | Para.ptrz = logical(Para.ptrz(NonVertPtr)); |
---|
887 | Para.Dist_Start = size(temp,1)+1; |
---|
888 | [x_ashIterator, alfa, status] = SigmoidLogBarrierSolver(Para, [], [], [], '', [], [], solverVerboseLevel); |
---|
889 | if any(S*x_ashIterator+inq > 0 ) |
---|
890 | disp('Inequality not satisfied'); |
---|
891 | max( S*x_ashIterator+inq) |
---|
892 | elseif status==2 |
---|
893 | fprintf([' Success with alfa=' num2str(alfa)]); |
---|
894 | end |
---|
895 | |
---|
896 | Para = x_ashIterator; |
---|
897 | |
---|
898 | tempPara = zeros(3*NuSubSupSize,1); |
---|
899 | tempPara(NonVertPtr) = Para; |
---|
900 | PlanePara(SubSupPtr) = tempPara; |
---|
901 | end |
---|
902 | end |
---|
903 | |
---|
904 | % porject the ray on planes to generate the FitDepth |
---|
905 | PlanePara = reshape(PlanePara,3,[]); |
---|
906 | FitDepth = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth); |
---|
907 | FitDepth(~maskSkyEroded) = (1./sum(PlanePara(:,Sup2Para(SupEpand(~maskSkyEroded))).*Ray(:,~maskSkyEroded),1))'; |
---|
908 | FitDepth = reshape(FitDepth,Default.VertYNuDepth,[]); |
---|
909 | % ==========Storage ============== |
---|
910 | if Default.Flag.AfterInferenceStorage |
---|
911 | save([ Default.ScratchFolder '/' strrep( Default.filename{1},'.jpg','') '_AInfnew.mat' ], 'FitDepth', 'depthMap', ... |
---|
912 | 'Sup', 'SupOri', 'RayOri','Ray','SupNeighborTable','maskSky','maskG','MultiScaleSupTable','WeiM','VarMapRaw'); |
---|
913 | end |
---|
914 | % =============================== |
---|
915 | [Position3DFited] = im_cr2w_cr(FitDepth,permute(Ray,[2 3 1])); |
---|
916 | Position3DFited(3,:) = -Position3DFited(3,:); |
---|
917 | Position3DFited = permute(Position3DFited,[2 3 1]); |
---|
918 | RR =permute(Ray,[2 3 1]); |
---|
919 | temp = RR(:,:,1:2)./repmat(RR(:,:,3),[1 1 2]); |
---|
920 | PositionTex = permute(temp./repmat(cat(3,Default.a_default,Default.b_default),... |
---|
921 | [Default.VertYNuDepth Default.HoriXNuDepth 1])... |
---|
922 | +repmat(cat(3,Default.Ox_default,Default.Oy_default),... |
---|
923 | [Default.VertYNuDepth Default.HoriXNuDepth 1]),[3 1 2]); |
---|
924 | PositionTex = permute(PositionTex,[2 3 1]); |
---|
925 | |
---|
926 | % write wrl file |
---|
927 | fprintf([' ' num2str( toc(inferenceTime) ) '\n : Writing WRL.']); |
---|
928 | WrlFacestHroiReduce(Position3DFited, PositionTex, SupOri, Default.filename{1}, Default.filename{1}, ... |
---|
929 | Default.OutPutFolder, GridFlag, 0); |
---|
930 | |
---|
931 | system(['gzip -9 -c ' Default.OutPutFolder Default.filename{1} '.wrl > ' ... |
---|
932 | Default.OutPutFolder Default.filename{1} '.wrl.gz']); |
---|
933 | copyfile([Default.OutPutFolder Default.filename{1} '.wrl.gz'], ... |
---|
934 | [Default.OutPutFolder Default.filename{1} '.wrl'],'f'); |
---|
935 | delete([Default.OutPutFolder Default.filename{1} '.wrl.gz']); |
---|
936 | |
---|
937 | return; |
---|