source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/LearningCode/Inference/OldVersion/ReportPlaneParaMRF_Conditioned_Ray.m @ 37

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

Added original make3d
File size: 57.1 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  ReportPlaneParaMRF_Conditioned(Default, step, DepthFolder,...
40          Sup, SupOri, MedSup,depthMap,VarMap,RayOri, Ray, SupNeighborTable,...
41          MedRay,maskSky,maskG,Algo,k, CornerList, OccluList,...
42          MultiScaleSupTable, StraightLineTable, HBrokeBook, VBrokeBook,previoslyStored,...
43          baseline);
44% This function runs the RMF over the plane parameter of each superpixels
45
46if nargin <20
47   baseline = 0;
48end
49solverVerboseLevel = 0;
50
51
52% ============== parameters for the decomposition problem
53% Optimal parameters for current code:
54% For one machine: Sedumi: (1,1)        Lpsolve:(3,1)
55% Multiple machines: Sedumi: (4,1)      LPsolve:(3,1)
56% lpsolve running time is 22 seconds for (4,2) arrangement; but numerical
57% accuracy needs to be resolved first.
58XNuDecompose = 3;       
59YNuDecompose = 1;
60% ============ parameters for the decomposition problem
61
62
63
64
65inferenceTime = tic;
66
67fprintf(['\n     : Building Matrices....       ']);
68
69% initialize parameters
70GridFlag = 0;
71NOYALMIP = 1;
72NoSecondStep = 0;
73Dual = false;
74displayFlag = false;
75RenderVrmlFlag = true;
76PlaneParaSegmentationFlag = false;
77
78scale = 1;
79StickHori = 5;  %0.1; % sticking power in horizontal direction
80StickVert = 10;     % sticking power in vertical direction
81Center = 2; % Co-Planar weight at the Center of each superpixel
82HoriConf = 1; % set the confidant of the learned depth at the middle in Horizontal direction of the image
83VertConf = 0.01; % set the confidant of the learned depth at the top of the image
84BandWith = 1; % Nov29 1 Nov30 0.1 check result change to 0.1 12/1 0.1 lost detail
85mapVert = linspace(VertConf,1,Default.VertYNuDepth); % modeling the gravity prior
86mapHori = [linspace(HoriConf,1,round(Default.HoriXNuDepth/2)) fliplr(linspace(HoriConf,1,Default.HoriXNuDepth-round(Default.HoriXNuDepth/2)))];
87
88% ========set the range of depth that our model in
89ClosestDist = 1;
90% set the FarestDist to very 5 times to median depth
91FarestDist = 1.5*median(depthMap(:)); % tried on university % nogood effect but keep it since usually it the rangelike this   % change to 1.5 for church
92% ================================================
93
94ceiling = 0*Default.VertYNuDepth; % set the position of the ceiling, related to No plane coming back constrain % changed for newchurch
95Name{1} = 'FraWOPri';
96Name{2} = 'FraCoP';
97if isempty(MultiScaleSupTable)
98   Name{3} = 'Var_FraStickCoP';
99else
100   Name{3} = 'Var_FraStickCoPSTasCoP';
101end
102if ~isempty(MultiScaleSupTable)
103   MultiScaleFlag = true;
104   WeiV = 2*ones(1,size(MultiScaleSupTable,2)-1);
105else
106   MultiScaleFlag = false;
107   WeiV = 1;
108end
109WeiV(1,1:2:end) = 6; % emphasize the middle scale three times smaller than large scale
110WeiV = WeiV./sum(WeiV);% normalize if pair of superpixels have same index in all the scale, their weight will be 10
111ShiftStick = -.1;  % between -1 and 0, more means more smoothing.
112ShiftCoP = -.5;  % between -1 and 0, more means more smoothing.
113gravity =true; % if true, apply the HoriConf and VertConf linear scale weight
114CoPST = true; % if true, apply the Straight line prior as the Co-Planar constrain
115
116%ConerImprove = false;
117%FractionalDepthError = true;
118
119
120% get rid of the error point and sky point in the depthMap
121% set every depth bigger than FarestDistmeter to FarestDistmeters
122%CleanedDepthMap = (depthMapif ~previoslyStored >80).*medfilt2(depthMap,[4 4])+(depthMap<=80).*depthMap;
123CleanedDepthMap = depthMap;
124%CleanedDepthMap(depthMap>FarestDist) = FarestDist; % don't clean the point >80 sometimes it occlusion
125%disp('Nu of depthMap>FarestDist');
126%sum(sum(depthMap>FarestDist))
127CleanedDepthMap(depthMap>FarestDist) = NaN; % don't clean the point >80 sometimes it occlusion
128Posi3D = im_cr2w_cr(CleanedDepthMap,permute(Ray,[2 3 1]));
129
130if ~previoslyStored
131
132   %NewMap = [rand(max(Sup(:)),3); [0 0 0]];
133   
134   % Clean the Sup near sky
135   maskSky = Sup == 0;
136   maskSkyEroded = imerode(maskSky, strel('disk', 4) );
137   SupEpand = ExpandSup2Sky(Sup,maskSkyEroded);
138   NuPatch = Default.HoriXNuDepth*Default.VertYNuDepth-sum(maskSky(:));
139
140   NuSup = setdiff(unique(Sup)',0);
141   NuSup = sort(NuSup);
142   NuSupSize = size(NuSup,2);
143
144   % Sup index and planeParameter index inverse map
145   Sup2Para = sparse(1,max(Sup(:)));
146   Sup2Para(NuSup) = 1:NuSupSize;
147
148   % constructinf the Straight line prior matrix Will be add in the CoPlane matrix
149   NuLine = size(StraightLineTable,2);
150   CoPSTList = [];
151
152   for i = 1:NuLine
153       L = StraightLineTable{i};
154       X = L(1:(end-1))';
155       Y = L(2:end)';
156       if isempty(X)
157          continue;
158       end
159       for j = 1:size(X,1)
160           if X(j)~=Y(j)
161              CoPSTList = [CoPSTList; X(j) Y(j)];
162           end
163       end
164   end
165end
166
167% Generate the Matrix for MRF
168 % ===========================================================================================================================================
169groundThreshold = cos([ zeros(1, Default.VertYNuDepth - ceil(Default.VertYNuDepth/2)+10) ...
170                                linspace(0,15,ceil(Default.VertYNuDepth/2)-10)]*pi/180);
171    %  v1 15 v2 20 too big v3 20 to ensure non misclassified as ground.
172    %  verticalThreshold = cos(linspace(5,55,Default.VertYNuDepth)*pi/180); % give a vector of size 55 in top to down :
173verticalThreshold = cos([ 5*ones(1,Default.VertYNuDepth - ceil(Default.VertYNuDepth/2)) ...
174                                linspace(5,55,ceil(Default.VertYNuDepth/2))]*pi/180);
175        % give a vector of size 55 in top to down :
176        % 50 means suface norm away from y axis more than 50 degree
177 % ===========================================================================================================================================
178
179PosiM = sparse(0,0);
180VarM = sparse(0,0);
181RayMd = sparse(0,0);
182RayAllOriM = sparse(0,0);
183RayAllM = sparse(0,0);
184RayMtilt = sparse(0,0);
185RayMCent = sparse(0,0);
186DepthInverseMCent = [];
187DepthInverseM = [];
188YPointer = [];
189YPosition = [];
190beta = [];
191EmptyIndex = [];
192
193for i = NuSup
194%    mask = Sup ==i;
195    mask = SupEpand ==i; % include the Ray that will be use to expand the NonSky
196    RayAllOriM = blkdiag( RayAllOriM, RayOri(:,mask)');
197    RayAllM = blkdiag( RayAllM, Ray(:,mask)');
198    mask = Sup ==i; % Not include the Ray that will be use to expand the NonSky   
199    [yt x] = find(mask);
200    CenterX = round(median(x));
201    CenterY = round(median(yt));
202    YPointer = [YPointer; CenterY >= ceiling];
203    YPosition = [YPosition; CenterY];
204    mask(isnan(CleanedDepthMap)) = false;
205    SupNuPatch(i) = sum(mask(:));
206%    if sum(mask(:)) < 5
207%       EmptyIndex = [EmptyIndex; i];
208%       mask(:) = false;
209%    end
210    % find center point
211    [yt x] = find(mask);
212    CenterX = round(median(x));
213    CenterY = round(median(yt));
214 
215  if ~all(mask(:)==0)
216    if gravity
217      if any(CleanedDepthMap(mask) <=0)
218         CleanedDepthMap(mask)
219%         pause
220      end
221      Pa2CenterRatio = median(CleanedDepthMap(mask))./CleanedDepthMap(mask);
222      if sum(mask(:)) > 0
223         RayMtilt = blkdiag(RayMtilt, ...
224             ( Pa2CenterRatio(:,[1 1 1]).*repmat(RayOri(:,CenterY,CenterX)',[ SupNuPatch(i) 1])- RayOri(:,mask)'));
225      else
226         RayMtilt = blkdiag(RayMtilt, RayOri(:,mask)');
227      end
228      RayMCent = blkdiag(RayMCent, RayOri(:,CenterY,CenterX)'*SupNuPatch(i)*mapVert(CenterY)*mapHori(CenterX));
229      PosiM = blkdiag(PosiM,Posi3D(:,mask)');%.*repmat( mapVert(yt)',[1 3]).*repmat( mapHori(x)',[1 3]));
230      VarM = [VarM; VarMap(mask).*(mapVert(yt)').*( mapHori(x)')];
231      RayMd = blkdiag(RayMd,RayOri(:,mask)'.*repmat( mapVert(yt)',[1 3]).*repmat( mapHori(x)',[1 3]));
232      DepthInverseMCent = [DepthInverseMCent; 1./median(CleanedDepthMap(mask)).*SupNuPatch(i)* mapVert(CenterY)'*mapHori(CenterX)];
233      DepthInverseM = [DepthInverseM; 1./CleanedDepthMap(mask).* mapVert(yt)'.*mapHori(x)'];
234    else
235      Pa2CenterRatio = CleanedDepthMap(CenterY,CenterX)./CleanedDepthMap(mask);
236      if sum(mask(:)) > 0
237         RayMtilt = blkdiag(RayMtilt, ...
238             ( Pa2CenterRatio(:,[1 1 1]).*repmat(RayOri(:,CenterY,CenterX)',[ SupNuPatch(i) 1])- RayOri(:,mask)'));
239      else
240         RayMtilt = blkdiag(RayMtilt, RayOri(:,mask)');
241      end
242      RayMCent = blkdiag(RayMCent, RayOri(:,CenterY,CenterX)'*SupNuPatch(i));
243      PosiM = blkdiag(PosiM,Posi3D(:,mask)');
244      VarM = [VarM; VarMap(mask)];
245      RayMd = blkdiag(RayMd,RayOri(:,mask)');
246      DepthInverseMCent = [DepthInverseMCent; 1./median(CleanedDepthMap(mask)).*SupNuPatch(i)];
247      DepthInverseM = [DepthInverseM; 1./CleanedDepthMap(mask)];
248    end
249  else
250     RayMtilt = blkdiag(RayMtilt, RayOri(:,mask)');
251     RayMCent = blkdiag(RayMCent, RayOri(:,mask)');
252     PosiM = blkdiag(PosiM, Posi3D(:,mask)');
253     VarM = [VarM; VarMap(mask)];
254     RayMd = blkdiag(RayMd, RayOri(:,mask)');
255  end
256end
257YPointer(YPointer==0) = -1;
258
259% buliding CoPlane Matrix=========================================================================
260
261NuSTList = 0;
262if CoPST
263   NuSTList = size(CoPSTList,1);
264   if ~isempty(CoPSTList)
265      [V H] = size(SupNeighborTable);
266      SupNeighborTable( CoPSTList(:,1)*V + CoPSTList(:,2)) = 1;
267      SupNeighborTable( CoPSTList(:,2)*V + CoPSTList(:,1)) = 1;
268%   ClosestNList = [ClosestNList; CoPSTList];
269   end
270end
271CoPM1 = sparse(0,3*NuSupSize);
272CoPM2 = sparse(0,3*NuSupSize);
273CoPEstDepth = sparse(0,0);
274CoPNorM = [];
275WeiCoP = [];
276
277for i = NuSup
278    Neighbor = find( SupNeighborTable(i,:) ~=0);
279    Neighbor = Neighbor( Neighbor> i);
280    for j = Neighbor
281        %  if ~CornerList(i)
282        mask = Sup == i;
283        SizeMaskAll = sum(mask(:));
284        [y x] = find(mask);
285        CenterX = round(median(x));
286        CenterY = round(median(y));
287        y = find(mask(:,CenterX));
288        if ~isempty(y)
289           CenterY = round(median(y));
290        end
291   
292        temp1 = sparse(1, 3*NuSupSize);
293        temp2 = sparse(1, 3*NuSupSize);
294        temp1(:,(Sup2Para( i)*3-2): Sup2Para( i)*3) = ...
295                                                   Ray(:,CenterY,CenterX)';
296        temp2(:,(Sup2Para( j)*3-2): Sup2Para( j)*3) = ...
297                                                   Ray(:,CenterY,CenterX)';
298    %NuNei-NuSTList
299
300%        if i < NuNei-NuSTList % only immediate connecting superpixels are weighted using MultiScaleSup
301%       wei = WeiV*10;%*(MultiScaleSupTable(Sup2Para(ClosestNList(i,1)),2:end) == MultiScaleSupTable(Sup2Para(ClosestNList(i,2)),2:end))'; 
302           if MultiScaleFlag
303              vector = (MultiScaleSupTable(Sup2Para( i),2:end) == MultiScaleSupTable(Sup2Para( j),2:end));
304              expV = exp(-10*(WeiV*vector' + ShiftCoP) );
305              wei = 1/(1+expV);
306           else
307              wei = 1;
308           end
309%        else
310%           wei = 1;
311%        end
312        oneRay1 = temp1*wei;
313        oneRay2 = temp2*wei;
314    %CoPM1 = [CoPM1; temp1*wei];
315    %CoPM2 = [CoPM2; temp2*wei];
316
317   
318        tempWeiCoP = [SizeMaskAll];
319        CoPEstDepth = [CoPEstDepth; max(median(CleanedDepthMap(mask)),ClosestDist)];
320   
321        mask = Sup == j;
322        SizeMaskAll = sum(mask(:));
323        [y x] = find(mask);
324        CenterX = round(median(x));
325        CenterY = round(median(y));
326        y = find(mask(:,CenterX));
327        if ~isempty(y)
328           CenterY = round(median(y));
329        end
330
331        temp1 = sparse(1, 3*NuSupSize);
332        temp2 = sparse(1, 3*NuSupSize);
333        temp1(:,(Sup2Para( i)*3-2): Sup2Para( i)*3) = ...
334                                                            Ray(:,CenterY,CenterX)';
335        temp2(:,(Sup2Para( j)*3-2): Sup2Para( j)*3) = ...
336                                                            Ray(:,CenterY,CenterX)';
337   
338    % Instead of having separate L-1 terms for symmetric co-planar constraint; do the following:
339    % If the penaly was ||.||_2^2 + ||.||_2^2; then the co-efficients are
340    % some kind of average of two rays.  For one norm; we take its average.
341    % (do not divide by 2 because the penalty should be double.
342   
343       CoPM1 = [CoPM1; temp1*wei + oneRay1 ];
344       CoPM2 = [CoPM2; temp2*wei + oneRay2 ];
345   
346       tempWeiCoP = [tempWeiCoP; SizeMaskAll];
347       WeiCoP = [WeiCoP; tempWeiCoP];
348       CoPEstDepth = [CoPEstDepth; max(median(CleanedDepthMap(mask)),ClosestDist)];
349  end
350end%=========================================================================================================
351
352%== find the boundary point that might need to be stick ot each other==========================================
353HoriStickM_i = sparse(0,3*NuSupSize);
354HoriStickM_j = sparse(0,3*NuSupSize);
355HoriStickPointInd = [];
356EstDepHoriStick = [];
357
358MAX_POINTS_STITCH_HORI = 2;    % the actual code will be modified in another copy of the file
359MIN_POINTS_STITCH = 2;      % ERROR:  not used.
360
361% ================================================================
362% NOTE: The actual algorithm should be picking precisely 2 points which are
363% FARTHEST away from the candidate set of neighbors.  This algorithm
364% will ALWAYS work and produce no surprising results.
365
366% In some cases, one may experiment with picking only 1 point when the 2
367% points are too close ---- this will make the algorithm faster; but might
368% produce surprising (e.g. a triangle sticking out) sometimes.
369% An ideal algorithm will reduce the number of points by checking for loops
370% passing through 3 or less superpixels through this matrix; and removing
371% them such that the smallest loop passes through 4 superpixels. (see EE263
372% for a quick algorithm to do this -- involves product of matrices.
373% ==================================================================
374
375DIST_STICHING_THRESHOLD_HORI = 0.4;   
376DIST_STICHING_THRESHOLD_HORI_ONLYCOL = -0.5;    % effectively not used,
377
378SupPixelNeighborList = sparse( max(Sup(:)), max(Sup(:)) );
379SupPixelParsedList = sparse( max(Sup(:)), max(Sup(:)) );
380recordAdded1 = sparse( max(Sup(:)), max(Sup(:)) );
381recordAdded2 = sparse( max(Sup(:)), max(Sup(:)) );
382addedIndexList = [ ];
383
384BounaryPHori = conv2(Sup,[1 -1],'same') ~=0;
385BounaryPHori(:,end) = 0;
386BounaryPVert = conv2(Sup,[1; -1],'same') ~=0;
387BounaryPVert(end,:) = 0;
388boundariesHoriIndex = find(BounaryPHori==1)';
389for i = boundariesHoriIndex
390    j = i+Default.VertYNuDepth;
391    if Sup(i) == 0 || Sup(j) == 0
392       continue;
393    end
394    SupPixelParsedList(Sup(i),Sup(j)) = SupPixelParsedList(Sup(i),Sup(j)) + 1;
395   
396    if SupPixelNeighborList(Sup(i),Sup(j)) == 0
397        recordAdded1(Sup(i),Sup(j)) = i;     
398    elseif SupPixelNeighborList(Sup(i),Sup(j)) >= MAX_POINTS_STITCH_HORI
399        continue;
400    elseif SupPixelNeighborList(Sup(i),Sup(j)) == 1  % inside this remove the close stiching terms
401        rowN = ceil(i/55);      colN = rem(i,55);
402        rowN_older = ceil( recordAdded1(Sup(i),Sup(j)) / 55);
403        colN_older = rem( recordAdded1(Sup(i),Sup(j)), 55);
404        if abs(rowN - rowN_older) + (55/305)*abs(colN - colN_older) > DIST_STICHING_THRESHOLD_HORI && ...
405            abs(colN - colN_older) > DIST_STICHING_THRESHOLD_HORI_ONLYCOL
406            recordAdded2(Sup(i),Sup(j)) = i;
407        else
408            continue;
409        end
410    elseif SupPixelNeighborList(Sup(i),Sup(j)) == 2     %Assuming MAX_POINTS_STITCH = 3
411        rowN = ceil(i/55);      colN = rem(i,55);
412        rowN_older1 = ceil( recordAdded1(Sup(i),Sup(j)) / 55);
413        colN_older1 = rem( recordAdded1(Sup(i),Sup(j)), 55);
414        rowN_older2 = ceil( recordAdded2(Sup(i),Sup(j)) / 55);
415        colN_older2 = rem( recordAdded2(Sup(i),Sup(j)), 55);
416       
417        if abs(rowN - rowN_older1) + (55/305)*abs(colN - colN_older1) > DIST_STICHING_THRESHOLD_HORI && ...
418            abs(rowN - rowN_older2) + (55/305)*abs(colN - colN_older2) > DIST_STICHING_THRESHOLD_HORI
419            ;
420        else
421            continue;
422        end   
423    end
424
425    % If you come here, it means you are probably adding it.
426    SupPixelNeighborList(Sup(i),Sup(j)) = SupPixelNeighborList(Sup(i),Sup(j)) + 1;
427    addedIndexList = [addedIndexList i];
428end
429
430
431WeightHoriNeighborStitch = [ ];
432
433for i = addedIndexList
434    j = i+Default.VertYNuDepth;
435   
436    WeightHoriNeighborStitch = [WeightHoriNeighborStitch;  SupPixelParsedList(Sup(i),Sup(j)) / ...
437                                            SupPixelNeighborList(Sup(i),Sup(j)) ];
438
439    Target(1) = Sup2Para(Sup(i));
440    Target(2) = Sup2Para(Sup(j));
441    rayBoundary(:,1) =  Ray(:,i);
442    rayBoundary(:,2) =  Ray(:,i);
443%    betaTemp = StickHori;%*(DistStickLengthNormWei.^2)*beta(Target(I));
444%    betaTemp = StickHori*WeiV;%*(MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end))';%*(DistStickLengthNormWei.^2)*beta(Target(I));
445    if MultiScaleFlag
446          vector = (MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end));
447          expV = exp(-10*(WeiV*vector' + ShiftStick) );
448       betaTemp = StickHori*(0.5+1/(1+expV)); %*(DistStickLengthNormWei.^2)*beta(Target(I));
449       % therr should always be sticking (know don't care about occlusion)
450    else
451       betaTemp = StickHori;
452    end
453    temp = sparse(3,NuSupSize);
454    temp(:,Target(1)) = rayBoundary(:,1);
455    HoriStickM_i = [HoriStickM_i; betaTemp*temp(:)'];
456    temp = sparse(3,NuSupSize);
457    temp(:,Target(2)) = rayBoundary(:,2);
458    HoriStickM_j = [HoriStickM_j; betaTemp*temp(:)'];
459    EstDepHoriStick = [EstDepHoriStick; sqrt(max(CleanedDepthMap(i),ClosestDist)*max(CleanedDepthMap(j),ClosestDist))];
460    HoriStickPointInd = [HoriStickPointInd i ];
461
462end
463% ==============================================
464
465
466
467% ======== finding the unucessary stiching points in Vertical direction ====
468VertStickM_i = sparse(0,3*NuSupSize);
469VertStickM_j = sparse(0,3*NuSupSize);
470VertStickPointInd = [];
471EstDepVertStick = [];
472
473MAX_POINTS_STITCH_VERT = 4;     %3
474DIST_STICHING_THRESHOLD_VERT = 0.1;     %0.3
475DIST_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.
476
477SupPixelNeighborList = sparse( max(Sup(:)), max(Sup(:)) );
478SupPixelParsedList = sparse( max(Sup(:)), max(Sup(:)) );
479recordAdded1 = sparse( max(Sup(:)), max(Sup(:)) );
480recordAdded2 = sparse( max(Sup(:)), max(Sup(:)) );
481addedIndexList = [ ];
482
483for i = find(BounaryPVert==1)'
484    j = i+1;
485    if Sup(i) == 0 || Sup(j) == 0
486       continue;
487    end
488    SupPixelParsedList(Sup(i),Sup(j)) = SupPixelParsedList(Sup(i),Sup(j)) + 1;
489   
490    if SupPixelNeighborList(Sup(i),Sup(j)) == 0
491        recordAdded1(Sup(i),Sup(j)) = i;
492    elseif SupPixelNeighborList(Sup(i),Sup(j)) >= MAX_POINTS_STITCH_VERT
493        continue;
494    elseif SupPixelNeighborList(Sup(i),Sup(j)) == 1  % inside this remove the close stiching terms
495        rowN = ceil(i/55);      colN = rem(i,55);
496        rowN_older = ceil( recordAdded1(Sup(i),Sup(j)) / 55);
497        colN_older = rem( recordAdded1(Sup(i),Sup(j)), 55);
498        if abs(rowN - rowN_older) + (55/305)*abs(colN - colN_older) > DIST_STICHING_THRESHOLD_VERT && ...
499                 abs(colN - colN_older) > DIST_STICHING_THRESHOLD_VERT_ONLYCOL
500            recordAdded2(Sup(i),Sup(j)) = i;
501        else
502            continue;
503        end
504    elseif SupPixelNeighborList(Sup(i),Sup(j)) == 2     %Assuming MAX_POINTS_STITCH = 3
505        rowN = ceil(i/55);      colN = rem(i,55);
506        rowN_older1 = ceil( recordAdded1(Sup(i),Sup(j)) / 55);
507        colN_older1 = rem( recordAdded1(Sup(i),Sup(j)), 55);
508        rowN_older2 = ceil( recordAdded2(Sup(i),Sup(j)) / 55);
509        colN_older2 = rem( recordAdded2(Sup(i),Sup(j)), 55);
510       
511        if abs(rowN - rowN_older1) + (55/305)*abs(colN - colN_older1) > DIST_STICHING_THRESHOLD_VERT && ...
512            abs(rowN - rowN_older2) + (55/305)*abs(colN - colN_older2) > DIST_STICHING_THRESHOLD_VERT
513            ;
514        else
515            continue;
516        end       
517     end
518
519    % If you come here, it means you are probably adding it.
520    SupPixelNeighborList(Sup(i),Sup(j)) = SupPixelNeighborList(Sup(i),Sup(j)) + 1;
521    addedIndexList = [addedIndexList i];
522end
523
524
525WeightVertNeighborStitch = [ ];
526
527for i = addedIndexList
528    j = i+1;
529   
530     WeightVertNeighborStitch = [WeightVertNeighborStitch;  SupPixelParsedList(Sup(i),Sup(j)) / ...
531                                            SupPixelNeighborList(Sup(i),Sup(j)) ];
532
533    Target(1) = Sup2Para(Sup(i));
534    Target(2) = Sup2Para(Sup(j));
535    rayBoundary(:,1) =  Ray(:,i);
536    rayBoundary(:,2) =  Ray(:,i);
537    if MultiScaleFlag
538       vector = (MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end));
539       expV = exp(-10*(WeiV*vector' + ShiftStick) );
540       betaTemp = StickVert*(0.5+1/(1+expV));
541       % therr should always be sticking (know don't care about occlusion)
542    else
543       betaTemp = StickVert;
544    end
545    temp = sparse(3,NuSupSize);
546    temp(:,Target(1)) = rayBoundary(:,1);
547    VertStickM_i = [VertStickM_i; betaTemp*temp(:)'];
548    temp = sparse(3,NuSupSize);
549    temp(:,Target(2)) = rayBoundary(:,2);
550    VertStickM_j = [VertStickM_j; betaTemp*temp(:)'];
551    EstDepVertStick = [EstDepVertStick; sqrt(max(CleanedDepthMap(i),ClosestDist)*max(CleanedDepthMap(j),ClosestDist))];
552    VertStickPointInd = [VertStickPointInd i ];
553
554end
555% ====finished finding the unucessary stiching points in Vertical direction ===
556
557
558% ======================================Finish building up matrix=====================hard work======================
559
560
561
562% ================================================================================================================
563depthfile = strrep(Default.filename{k},'img','depth_learned'); %
564
565
566% Start Decompose the image align with superpixels ==================================================================================
567% define the decomposition in X direction only
568
569
570 
571TotalRectX = 2*XNuDecompose-1;
572TotalRectY= 2*YNuDecompose-1;
573PlanePara = NaN*ones(3*NuSupSize,1); % setup the lookuptable for the solved plane parameter
574
575 opt = sdpsettings('solver','sedumi','cachesolvers',1,'sedumi.eps',1e-6,'verbose',solverVerboseLevel);
576%     opt = sdpsettings('solver','lpsolve','cachesolvers',1,'verbose',5);
577%     opt = sdpsettings('solver','lpsolve','cachesolvers',1,'showprogress',1, 'verbose',4);
578%    opt = sdpsettings('solver','glpk','cachesolvers',1, 'usex0', 1);   
579%    opt = sdpsettings('solver','glpk', 'usex0', 1);   
580    %opt = sdpsettings('solver','linprog','cachesolvers',1);   
581%    opt = sdpsettings('solver','bpmpd');
582   
583for k = 0:(TotalRectX-1)
584  l = rem(k*2,(TotalRectX));
585  RangeX = (1+ceil(Default.HoriXNuDepth/XNuDecompose)*l/2):...
586            min((1+ceil(Default.HoriXNuDepth/XNuDecompose)*(l/2+1)),Default.HoriXNuDepth);
587  RangeX = ceil(RangeX);
588  for q = 0:(TotalRectY-1)
589    l = rem(q*2,(TotalRectY));
590    RangeY = (1+ceil(Default.VertYNuDepth/YNuDecompose)*l/2):...
591              min((1+ceil(Default.VertYNuDepth/YNuDecompose)*(l/2+1)),Default.VertYNuDepth);
592    RangeY = ceil(RangeY);
593    mask = zeros(size(Sup));
594    mask(RangeY,RangeX) = 1;
595    mask =logical(mask);
596    SubSup = sort(setdiff(unique( reshape( Sup(RangeY,RangeX),1,[])),0));
597    BoundarySup = [];
598%    for m = SubSup
599    if any(SubSup <=0)
600       SubSup(SubSup<=0)
601    end
602    BoundarySup = find(sum(SupNeighborTable(SubSup,:), 1) ~=0);
603%        maskList = ClosestNList(:,1) == m;
604%        BoundarySup = [ BoundarySup ClosestNList(maskList,2)'];
605%    end
606    BoundarySup = unique(setdiff(BoundarySup,[0 SubSup] ));
607    % chech if BoundarySup non-NaN in PlanePara
608    checkNoNNaN = ~isnan(PlanePara(Sup2Para(BoundarySup)*3));
609    BoundarySup = BoundarySup(checkNoNNaN);
610    TotalSup = sort([SubSup BoundarySup]);
611
612    SubSupPtr = [ Sup2Para(SubSup)*3-2;...
613                    Sup2Para(SubSup)*3-1;...
614                    Sup2Para(SubSup)*3];
615    SubSupPtr = SubSupPtr(:);
616%    TotalSupPtr = [ Sup2Para(TotalSup)*3-2;...
617%                    Sup2Para(TotalSup)*3-1;...
618%                    Sup2Para(TotalSup)*3];
619%    TotalSupPtr = TotalSupPtr(:);
620    BoundarySupPtr = [ Sup2Para(BoundarySup)*3-2;...
621                       Sup2Para(BoundarySup)*3-1;...
622                       Sup2Para(BoundarySup)*3];
623    BoundarySupPtr = BoundarySupPtr(:);
624%    NuSubSupSize = size(TotalSup,2);
625    NuSubSupSize = size(SubSup,2);
626%    TotalSup2Para = sparse(1,max(TotalSup));
627%    TotalSup2Para(TotalSup) = 1:NuSubSupSize;
628%    BoundarySupPtrSub = [ TotalSup2Para(BoundarySup)*3-2;...
629%                          TotalSup2Para(BoundarySup)*3-1;...
630%                          TotalSup2Para(BoundarySup)*3];
631%    BoundarySupPtrSub = BoundarySupPtrSub(:);
632   
633    % clearn RayAllM PosiM CoPM1 HoriStickM_i VertStickM_i
634%    NewRayAllM = RayAllM(:,TotalSupPtr);
635    NewRayAllM = RayAllM(:,SubSupPtr);
636    tar = sum(NewRayAllM ~= 0,2) == 3;
637    NewRayAllM = NewRayAllM(tar,:);
638%    NewPosiM = PosiM(:,TotalSupPtr);
639    NewPosiM = PosiM(:,SubSupPtr);
640    tar = sum(NewPosiM ~= 0,2) == 3;
641    NewPosiM = NewPosiM(tar,:);
642    NewVarM = VarM(tar);
643   
644%    NewCoPM = CoPM1(:,TotalSupPtr) - CoPM2(:,TotalSupPtr);
645    NewCoPM = CoPM1(:,SubSupPtr) - CoPM2(:,SubSupPtr);
646    NewCoPMBound = CoPM1(:,BoundarySupPtr) - CoPM2(:,BoundarySupPtr);
647%    tar = sum(NewCoPM ~= 0,2) ==6;
648    tar = sum(NewCoPM ~= 0,2) + sum(NewCoPMBound ~= 0,2)==6;
649    NewCoPMBound = NewCoPMBound*PlanePara(BoundarySupPtr);
650    NewCoPM = NewCoPM(tar,:);
651    NewCoPMBound = NewCoPMBound(tar);
652    NewCoPEstDepth = CoPEstDepth(tar);   
653
654%    NewHoriStickM = HoriStickM_i(:,TotalSupPtr)-HoriStickM_j(:,TotalSupPtr);
655    NewHoriStickM = HoriStickM_i(:,SubSupPtr)-HoriStickM_j(:,SubSupPtr);
656    NewHoriStickMBound = HoriStickM_i(:,BoundarySupPtr)-HoriStickM_j(:,BoundarySupPtr);
657%    tar = sum(NewHoriStickM ~= 0,2) ==6;
658    tar = sum(NewHoriStickM ~= 0,2)+  sum(NewHoriStickMBound ~= 0,2) ==6;
659    NewHoriStickM = NewHoriStickM(tar,:);   
660%    TotalSupPtr = [ Sup2Para(TotalSup)*3-2;...
661%                    Sup2Para(TotalSup)*3-1;...
662%                    Sup2Para(TotalSup)*3];
663%    TotalSupPtr = TotalSupPtr(:);
664    BoundarySupPtr = [ Sup2Para(BoundarySup)*3-2;...
665                       Sup2Para(BoundarySup)*3-1;...
666                       Sup2Para(BoundarySup)*3];
667    BoundarySupPtr = BoundarySupPtr(:);
668%    NuSubSupSize = size(TotalSup,2);
669    NuSubSupSize = size(SubSup,2);
670%    TotalSup2Para = sparse(1,max(TotalSup));
671%    TotalSup2Para(TotalSup) = 1:NuSubSupSize;
672%    BoundarySupPtrSub = [ TotalSup2Para(BoundarySup)*3-2;...
673%                          TotalSup2Para(BoundarySup)*3-1;...
674%                          TotalSup2Para(BoundarySup)*3];
675%    BoundarySupPtrSub = BoundarySupPtrSub(:);
676   
677    % clearn RayAllM PosiM CoPM1 HoriStickM_i VertStickM_i
678%    NewRayAllM = RayAllM(:,TotalSupPtr);
679    NewRayAllM = RayAllM(:,SubSupPtr);
680    tar = sum(NewRayAllM ~= 0,2) == 3;
681    NewRayAllM = NewRayAllM(tar,:);
682%    NewPosiM = PosiM(:,TotalSupPtr);
683    NewPosiM = PosiM(:,SubSupPtr);
684    tar = sum(NewPosiM ~= 0,2) == 3;
685    NewPosiM = NewPosiM(tar,:);
686    NewVarM = VarM(tar);
687   
688%    NewCoPM = CoPM1(:,TotalSupPtr) - CoPM2(:,TotalSupPtr);
689    NewCoPM = CoPM1(:,SubSupPtr) - CoPM2(:,SubSupPtr);
690    NewCoPMBound = CoPM1(:,BoundarySupPtr) - CoPM2(:,BoundarySupPtr);
691%    tar = sum(NewCoPM ~= 0,2) ==6;
692    tar = sum(NewCoPM ~= 0,2) + sum(NewCoPMBound ~= 0,2)==6;
693    NewCoPMBound = NewCoPMBound*PlanePara(BoundarySupPtr);
694    NewCoPM = NewCoPM(tar,:);
695    NewCoPMBound = NewCoPMBound(tar);
696    NewCoPEstDepth = CoPEstDepth(tar);   
697
698%    NewHoriStickM = HoriStickM_i(:,TotalSupPtr)-HoriStickM_j(:,TotalSupPtr);
699    NewHoriStickM = HoriStickM_i(:,SubSupPtr)-HoriStickM_j(:,SubSupPtr);
700    NewHoriStickMBound = HoriStickM_i(:,BoundarySupPtr)-HoriStickM_j(:,BoundarySupPtr);
701%    tar = sum(NewHoriStickM ~= 0,2) ==6;
702    tar = sum(NewHoriStickM ~= 0,2)+  sum(NewHoriStickMBound ~= 0,2) ==6;
703    NewHoriStickM = NewHoriStickM(tar,:);   
704    NewEstDepHoriStick = EstDepHoriStick(tar);
705    NewHoriStickMBound = NewHoriStickMBound*PlanePara(BoundarySupPtr);
706    NewHoriStickMBound = NewHoriStickMBound(tar);
707        NewWeightHoriNeighborStitch = WeightHoriNeighborStitch(tar);
708
709%    NewVertStickM = VertStickM_i(:,TotalSupPtr)-VertStickM_j(:,TotalSupPtr);
710    NewVertStickM = VertStickM_i(:,SubSupPtr)-VertStickM_j(:,SubSupPtr);
711    NewVertStickMBound = VertStickM_i(:,BoundarySupPtr)-VertStickM_j(:,BoundarySupPtr);
712%    tar = sum(NewVertStickM ~= 0,2) ==6;
713    tar = sum(NewVertStickM ~= 0,2) + sum(NewVertStickMBound ~= 0,2)==6;
714    NewVertStickM = NewVertStickM(tar,:);
715    NewEstDepVertStick = EstDepVertStick(tar);
716    NewVertStickMBound = NewVertStickMBound*PlanePara(BoundarySupPtr);
717    NewVertStickMBound = NewVertStickMBound(tar);
718        NewWeightVertNeighborStitch = WeightVertNeighborStitch(tar);
719   
720 
721    for i = step
722        ParaPPCP = sdpvar(3*NuSubSupSize,1);
723        F = set(ParaPPCP(3*(1:NuSubSupSize)-1).*YPointer(Sup2Para(SubSup))<=0)...
724            +set(NewRayAllM*ParaPPCP <=1/ClosestDist)...
725            +set(NewRayAllM*ParaPPCP >=1/FarestDist);
726%            +set(ParaPPCP(BoundarySupPtrSub) == PlanePara(BoundarySupPtr) ); % Special constrain for decomp fix the solved neighbor plane parameter       
727%         +set(RayAllOriM(:,TotalSupPtr)*ParaPPCP >=1/FarestDist)...
728%         +set(RayAllOriM(:,TotalSupPtr)*ParaPPCP <=1/ClosestDist);
729% First fit the plane to find the estimated plane parameters
730% If PosiM contain NaN data the corresponding Plane Parameter will also be NaN
731    if i == 1 % W/O  priors
732%    solvesdp(F,norm( (PosiM*ParaPPCP-ones(size(PosiM,1),1))./(abs(VarM)),1)...
733       solvesdp(F,norm( (PosiM*ParaPPCP-ones(size(PosiM,1),1))./exp(abs(VarM)/BandWith),1)...
734             , opt);
735    elseif i ==2 % Coner
736       solvesdp(F,norm( (PosiM*ParaPPCP-ones(size(PosiM,1),1))./exp(abs(VarM)/BandWith),1)...
737             +Center*norm(((CoPM1 - CoPM2)*ParaPPCP).*CoPEstDepth, 1)...
738             , opt);
739    else % all the priors
740%    sum(sum(isnan(PosiM)))
741       
742%     solvesdp(F,norm( (PosiM*ParaPPCP-ones(size(PosiM,1),1))./exp(abs(VarM)/BandWith),1)...
743%              +Center*norm(((CoPM1 - CoPM2)*ParaPPCP).*CoPEstDepth, 1)...
744%              +norm(((HoriStickM_i-HoriStickM_j)*ParaPPCP).*EstDepHoriStick,1)...
745%              +norm(((VertStickM_i-VertStickM_j)*ParaPPCP).*EstDepVertStick,1)...
746%              , opt);
747   
748
749% =========  Debugging the code for numerical accuracy  ==========
750
751
752
753% ================================================================
754
755        WeightsSelfTerm = 1 ./ exp(abs(NewVarM)/BandWith);
756   
757   
758    %============ solver specific options ======
759    if strcmp(opt.solver, 'glpk')
760        disp('Using x0 for GPLK');
761        ParaPPCP = NewPosiM \ ones(size(NewPosiM,1),1);
762    end
763   fprintf(['     ' num2str( toc(inferenceTime) ) '\n     : In 1st level Optimization, using ' opt.solver '.  ' ...
764            '(' num2str(k+1) '/' num2str((TotalRectX-1)+1) ',' num2str(l+1) '/' num2str((TotalRectY-1)+1) ')']);
765
766       if NOYALMIP
767          sol = solvesdp(F, norm( (NewPosiM*ParaPPCP-ones(size(NewPosiM,1),1)) .* WeightsSelfTerm,1)...
768                        + Center*norm(((NewCoPM)*ParaPPCP + NewCoPMBound).*NewCoPEstDepth, 1)...
769                        + norm(((NewHoriStickM)*ParaPPCP + NewHoriStickMBound).*...
770                                 NewEstDepHoriStick.*NewWeightHoriNeighborStitch,1)...
771                        + norm(((NewVertStickM)*ParaPPCP + NewVertStickMBound).*...
772                                 NewEstDepVertStick.*NewWeightVertNeighborStitch,1) ...
773                        , opt);
774%opt
775%mean(abs(ParaPPCP(:))) 
776
777     else
778       
779       
780
781       % ============Using Sedumi directly ==================================
782          if Dual
783          % trial for dual form
784
785
786          else
787            % trial for MS_E212 method: Primal form
788            %  D1 = PosiM;
789            %  D2 = (CoPM1 - CoPM2);
790            %  D3 = (HoriStickM_i-HoriStickM_j);
791            %  D4 = (VertStickM_i-VertStickM_j);
792            %  b = -ones(size(PosiM,1),1);
793 
794            %  C1 = RayAllM;
795            %  C2 = -RayAllM
796            %  clear tempSparse;
797            %  e1 = 1/ClosestDist*ones(size(RayAllM,1),1);
798            %  e2 = -1/ FarestDist*ones(size(RayAllM,1),1);
799
800
801            %  e2 = sparse(3*NuSupSize,1);
802            %  C2 = spdiags(tempSparse,0,3*NuSupSize,3*NuSupSize);
803
804            %  t1 = 1./exp(abs(VarM)/BandWith);
805            %  t2 = CoPEstDepth*Center;
806            %  t3 = EstDepHoriStick;
807            %  t4 = EstDepVertStick;
808
809            sizeX = size(PosiM,2);
810            sizeD1 = size(PosiM,1);
811            sizeD2 = size(CoPM1,1);
812            sizeD3 = size(HoriStickM_i,1);
813            sizeD4 = size(VertStickM_i,1);
814            sizeC1 = size(RayAllM,1);
815            sizeC2 = sizeC1;
816            %   sizeC2 = size(spdiags(tempSparse,0,3*NuSupSize,3*NuSupSize),1);
817 
818 
819            A = [ -2*PosiM ... % x
820                  -speye(sizeD1,sizeD1) ... % V_1
821                  sparse(sizeD1,(sizeD2+sizeD3+sizeD4) );...
822                  ...
823                  -2*(CoPM1 - CoPM2) ... % x
824                  sparse(sizeD2,sizeD1) ...
825                  -spdiags(ones(sizeD2,1),0,sizeD2,sizeD2) ... % V_2
826                  sparse(sizeD2, (sizeD3+sizeD4));...
827                  ...
828                  -2*(HoriStickM_i-HoriStickM_j) ... %x
829                  sparse(sizeD3, (sizeD1+sizeD2) ) ...
830                  -spdiags(ones(sizeD3,1),0,sizeD3,sizeD3) ... % V_3
831                  sparse(sizeD3,sizeD4 );...
832                  ...
833                  -2*(VertStickM_i-VertStickM_j) ... %x
834                  sparse(sizeD4, (sizeD1+sizeD2+sizeD3)) ...
835                  -spdiags(ones(sizeD4,1),0,sizeD4,sizeD4); ... % V_4
836                  ...
837                  RayAllM ...
838                  sparse(sizeC1,(sizeD1+sizeD2+sizeD3+sizeD4) ); ...
839                  ...
840                  -RayAllM ...
841                  sparse(sizeC2,(sizeD1+sizeD2+sizeD3+sizeD4) ); ...
842            ]; 
843
844            b = [ -2*ones(size(PosiM,1),1); ...
845                  sparse(sizeD2+sizeD3+sizeD4,1);...
846                  1/ClosestDist*ones(size(RayAllM,1),1);...
847                  -1/FarestDist*ones(size(RayAllM,1),1);...
848                ];
849            H = [ (1./exp(abs(VarM)/BandWith))'*(PosiM)+...
850                  (CoPEstDepth*Center)'*(CoPM1 - CoPM2)+...
851                  (EstDepHoriStick)'*((HoriStickM_i-HoriStickM_j))+...
852                  (EstDepVertStick)'*((VertStickM_i-VertStickM_j)) ...
853                  ...
854                  (1./exp(abs(VarM)/BandWith))' ...
855                  ...
856                  CoPEstDepth'*Center ...
857                  ...
858                  EstDepHoriStick' ...
859                  ...
860                  EstDepVertStick' ...
861                ];
862            lb = [-Inf*ones(sizeX,1); ...
863                  sparse((sizeD1+sizeD2+sizeD3+sizeD4),1)];
864            tempSparse = sparse(3*NuSupSize,1);
865            tempSparse(3*(1:NuSupSize)-1) = YPointer;
866            ub = Inf*ones(sizeX,1);
867            ub(logical(tempSparse)) = 0;
868            ub = [ub; Inf*ones((sizeD1+sizeD2+sizeD3+sizeD4),1)];
869 
870            %   options = optimset('LargeScale','on');
871            %  K.f = sizeX;
872            %  ParaPPCP = J*sedumi(A,b,H,K);
873            %   linOut = linprog(H',A,b,[],[],lb,ub,[],options);
874            [obj,x,duals,stat] = lp_solve(-H',A,b,-1*ones(size(b)),lb,ub);
875             ParaPPCP = x(1:sizeX);
876            %   opt = sdpsettings('solver','sedumi');
877            %   ParaPPCP = sdpvar(sizeX + sizeD1 +sizeD2+sizeD3+sizeD4);
878            %   F = set(A*ParaPPCP <= b) + set(ParaPPCP>=lb);
879            %   ParaPPCP = J*solvesdp(F,H*ParaPPCP,opt);
880          end
881       end
882    end
883
884    ParaPPCP = double(ParaPPCP);
885%     sum(isnan(ParaPPCP))
886    yalmip('clear');
887    PlanePara(SubSupPtr) = ParaPPCP;
888
889% %    SepPointMeasureHori = 1./(HoriStickM_i*ParaPPCP)-1./(HoriStickM_j*ParaPPCP);
890% %    SepPointMeasureVert = 1./(VertStickM_i*ParaPPCP)-1./(VertStickM_j*ParaPPCP);
891%     ParaPPCP = reshape(ParaPPCP,3,[]);
892%     %any(any(isnan(ParaPPCP)))
893%     % porject the ray on planes to generate the ProjDepth
894%     FitDepthPPCP = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth);
895%     FitDepthPPCP(~maskSkyEroded & mask) = (1./sum(ParaPPCP(:,TotalSup2Para(SupEpand(~maskSkyEroded & mask))).*Ray(:,~maskSkyEroded & mask),1))';
896%     FitDepthPPCP = reshape(FitDepthPPCP,Default.VertYNuDepth,[]);
897%     % storage for comparison purpose ======================
898% %    depthMap = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth);
899% %    depthMap(~maskSkyEroded | mask) = (1./sum(ParaPPCP(:,TotalSup2Para(SupEpand(~maskSkyEroded | mask))).*RayOri(:,~maskSkyEroded | mask),1))';
900% %    depthMap = reshape(depthMap,Default.VertYNuDepth,[]);
901% %    if baseline == 0
902% %       system(['mkdir ' Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '/']);
903% %       save([Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '/' depthfile '.mat'],'depthMap');
904% %    elseif baseline == 1
905% %       system(['mkdir ' Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline/']);
906% %       save([Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline/' depthfile '.mat'],'depthMap');
907% %    else
908% %       system(['mkdir ' Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline2/']);
909% %       save([Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline2/' depthfile '.mat'],'depthMap');
910% %    end
911%   
912%     % =====================================================
913%     [Position3DFitedPPCP] = im_cr2w_cr(FitDepthPPCP,permute(Ray,[2 3 1]));
914%     if true %RenderVrmlFlag && i==3
915%          Position3DFitedPPCP(3,:) = -Position3DFitedPPCP(3,:);
916%          Position3DFitedPPCP = permute(Position3DFitedPPCP,[2 3 1]);
917%          RR =permute(Ray,[2 3 1]);
918%          temp = RR(:,:,1:2)./repmat(RR(:,:,3),[1 1 2]);
919%          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]);
920%          PositionTex = permute(PositionTex,[2 3 1]);
921%          WrlFacest(Position3DFitedPPCP,PositionTex,Sup,Default.filename{1},'./')
922%
923% %        [VrmlName] = vrml_test_faceset_goodSkyBoundary(      Default.filename{k}, Position3DFitedPPCP, FitDepthPPCP, permute(Ray,[2 3 1]), [Name{i} DepthFolder 'ExpVar'], ...
924% %                     [], [], 0, maskSkyEroded, maskG, 1, 0, Default.a_default, Default.b_default, Default.Ox_default, Default.Oy_default);
925% %        system(['gzip -9 -c ' Default.ScratchDataFolder '/vrml/' VrmlName ' > ' Default.ScratchDataFolder '/vrml/' VrmlName '.gz']);
926%        %delete([Default.ScratchDataFolder '/vrml/' VrmlName]);
927%     end
928%    save([Default.ScratchDataFolder '/data/PreDepth/FitDepthPPCP' num2str(k) '.mat'],'FitDepthPPCP','SepPointMeasureHori',...
929%      'SepPointMeasureVert','VertStickPointInd','HoriStickPointInd');
930%save([Default.ScratchDataFolder '/data/All.mat']);
931  end
932  end
933end
934% build the whole image
935    PlanePara = reshape(PlanePara,3,[]);
936    %any(any(isnan(PlanePara)))
937    % porject the ray on planes to generate the ProjDepth
938    FitDepthPPCP = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth);
939    FitDepthPPCP(~maskSkyEroded) = (1./sum(PlanePara(:,Sup2Para(SupEpand(~maskSkyEroded ))).*Ray(:,~maskSkyEroded ),1))';
940    FitDepthPPCP = reshape(FitDepthPPCP,Default.VertYNuDepth,[]);
941    % storage for comparison purpose ======================
942%    depthMap = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth);
943%    depthMap(~maskSkyEroded | mask) = (1./sum(PlanePara(:,TotalSup2Para(SupEpand(~maskSkyEroded | mask))).*RayOri(:,~maskSkyEroded | mask),1))';
944%    depthMap = reshape(depthMap,Default.VertYNuDepth,[]);
945%    if baseline == 0
946%       system(['mkdir ' Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '/']);
947
948
949
950%       save([Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '/' depthfile '.mat'],'depthMap');
951%    elseif baseline == 1
952%       system(['mkdir ' Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline/']);
953%       save([Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline/' depthfile '.mat'],'depthMap');
954%    else
955%       system(['mkdir ' Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline2/']);
956%       save([Default.ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline2/' depthfile '.mat'],'depthMap');
957%    end
958
959    % =====================================================
960    [Position3DFitedPPCP] = im_cr2w_cr(FitDepthPPCP,permute(Ray,[2 3 1]));
961    if NoSecondStep %RenderVrmlFlag && i==3
962         Position3DFitedPPCP(3,:) = -Position3DFitedPPCP(3,:);
963         Position3DFitedPPCP = permute(Position3DFitedPPCP,[2 3 1]);
964         RR =permute(Ray,[2 3 1]);
965         temp = RR(:,:,1:2)./repmat(RR(:,:,3),[1 1 2]);
966         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]);
967         PositionTex = permute(PositionTex,[2 3 1]);
968         WrlFacestHroiReduce(Position3DFitedPPCP,PositionTex,SupOri, [ Default.filename{1} '1st'],[ Default.filename{1} '1st'], ...
969                             Default.OutPutFolder, GridFlag, 0);
970         system(['gzip -9 -c ' Default.OutPutFolder Default.filename{1} '1st.wrl > ' ...
971                Default.OutPutFolder Default.filename{1} '1st.wrl.gz']);
972         system(['cp '  Default.OutPutFolder Default.filename{1} '1st.wrl.gz ' ...
973                Default.OutPutFolder Default.filename{1} '1st.wrl']);
974         delete([Default.OutPutFolder Default.filename{1} '1st.wrl.gz']);
975    end
976
977    % Trial for segmentaion of the plane parameters
978    if PlaneParaSegmentationFlag
979       TempSup = Sup;%imresize(Sup,[330 610]);
980       TempSup(TempSup==0) = max(TempSup(:))+1;
981       TempSup2Para = Sup2Para;
982       TempSup2Para(end+1) = size(PlanePara,2);
983       TempPlanePara = PlanePara;
984       TempPlanePara(:,end+1) = 0;
985       PlaneParaPics = TempPlanePara(:,TempSup2Para( TempSup));
986       PlaneParaPics = PlaneParaPics./repmat( 2*max(abs(PlaneParaPics),[],2), 1, size(PlaneParaPics,2));
987       PlaneParaPics = permute( reshape( PlaneParaPics, 3, 55, []), [2 3 1])+0.5;
988       MergedPlaneParaPics = segmentImgOpt( Default.sigm*scale, Default.k*scale, Default.minp*10*scale, PlaneParaPics, '', 0) + 1;
989       if Default.Flag.DisplayFlag
990          figure(400);
991          imaegsc( PlaneParaPics);
992          imaegsc( MergedPlaneParaPics);
993       end
994       save('/afs/cs/group/reconstruction3d/scratch/testE/PlaneParaPics.mat','PlaneParaPics','MergedPlaneParaPics');
995    end
996
997%return;
998%==================Finished for one step MRF==========================================================================================================
999
1000if NoSecondStep
1001        return;
1002end
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012% ====================================following are 2nd step MRF to give more visually pleasing result=======================================
1013;% generating new PosiMPPCP using the new losition
1014
1015% save([Default.ScratchDataFolder '/data/VertShiftVrml.mat' ] );
1016%  groundThreshold = cos(5*pi/180);  % make small range
1017%  verticalThreshold = cos(50*pi/180);
1018  normPara = norms(PlanePara);
1019  normalizedPara = PlanePara ./ repmat( normPara, [3 1]);
1020  groundPara = abs(normalizedPara(2,:)) >= groundThreshold(YPosition);
1021  groundParaInd = find(groundPara);
1022  verticalPara = abs(normalizedPara(2,:)) <= verticalThreshold(YPosition); % change to have different range of vertical thre ================
1023  verticalParaInd = find(verticalPara);
1024
1025  indexVertical = find( verticalPara)*3-1;
1026  indexGroundX = find( groundPara)*3-2;
1027  indexGroundZ = find( groundPara)*3;
1028
1029  PosiMPPCP = sparse(0,0);
1030  VarM2 = sparse(0,0);
1031  %  VertVar = 10^(-2);
1032  % forming new supporting matrix using new depth and get rid of the support of the vertical plane
1033  for i = NuSup
1034      mask = Sup == i;
1035      if any(verticalParaInd == Sup2Para(i))
1036         mask = logical(zeros(size(Sup)));
1037      end
1038         PosiMPPCP = blkdiag(PosiMPPCP, Posi3D(:,mask)');
1039         VarM2 = [VarM2; VarMap(mask)];         
1040  end
1041
1042% Start Decompose image =========================
1043  TotalRectX = 2*XNuDecompose-1;
1044  TotalRectY = 2*YNuDecompose-1;
1045  PlanePara = NaN*ones(3*NuSupSize,1); % setup the lookuptable for the solved plane parameter
1046
1047  opt = sdpsettings('solver','sedumi','cachesolvers',1, 'verbose', solverVerboseLevel);
1048  %     opt = sdpsettings('solver','lpsolve','cachesolvers',1);
1049  %    opt = sdpsettings('solver','glpk','cachesolvers',1);   
1050
1051  for k = 0:(TotalRectX-1)
1052    l = rem(k*2,(TotalRectX));
1053    RangeX = (1+ceil(Default.HoriXNuDepth/XNuDecompose)*l/2):...
1054              min((1+ceil(Default.HoriXNuDepth/XNuDecompose)*(l/2+1)),Default.HoriXNuDepth);
1055    RangeX = ceil(RangeX);
1056    for q = 0:(TotalRectY-1)
1057      l = rem(q*2,(TotalRectY));
1058      RangeY = (1+ceil(Default.VertYNuDepth/YNuDecompose)*l/2):...
1059                min((1+ceil(Default.VertYNuDepth/YNuDecompose)*(l/2+1)),Default.VertYNuDepth);
1060      RangeY = ceil(RangeY);
1061      mask = zeros(size(Sup));
1062      mask(RangeY,RangeX) = 1;
1063      mask =logical(mask);
1064%       SubSup = sort(setdiff(reshape( Sup(RangeY,RangeX),[],1),0))';
1065      SubSup = sort(setdiff( unique( reshape( Sup(RangeY,RangeX),1,[])),0));
1066      BoundarySup = [];
1067%      for m = SubSup
1068      BoundarySup = find(sum(SupNeighborTable(SubSup,:), 1) ~=0);
1069%          maskList = ClosestNList(:,1) == m;
1070%          BoundarySup = [ BoundarySup ClosestNList(maskList,2)'];
1071%      end
1072      BoundarySup = unique(setdiff(BoundarySup,[0 SubSup] ));
1073
1074      % chech if BoundarySup non-NaN in PlanePara
1075      checkNoNNaN = ~isnan(PlanePara(Sup2Para(BoundarySup)*3));
1076      BoundarySup = BoundarySup(checkNoNNaN);
1077      TotalSup = sort([SubSup BoundarySup]);
1078
1079%      TotalSupPtr = [ Sup2Para(TotalSup)*3-2;...
1080%                      Sup2Para(TotalSup)*3-1;...
1081%                      Sup2Para(TotalSup)*3];
1082%      TotalSupPtr = TotalSupPtr(:);
1083      SubSupPtr = [ Sup2Para(SubSup)*3-2;...
1084                    Sup2Para(SubSup)*3-1;...
1085                    Sup2Para(SubSup)*3];
1086      SubSupPtr = SubSupPtr(:);
1087      BoundarySupPtr = [ Sup2Para(BoundarySup)*3-2;...
1088                         Sup2Para(BoundarySup)*3-1;...
1089                         Sup2Para(BoundarySup)*3];
1090      BoundarySupPtr =BoundarySupPtr(:);
1091%      NuSubSupSize = size(TotalSup,2);
1092      NuSubSupSize = size(SubSup,2);
1093%      TotalSup2Para = sparse(1,max(TotalSup));
1094%      TotalSup2Para(TotalSup) = 1:NuSubSupSize;
1095      SubSup2Para = sparse(1,max(SubSup));
1096      SubSup2Para(SubSup) = 1:NuSubSupSize;
1097%      BoundarySupPtrSub = [ TotalSup2Para(BoundarySup)*3-2;...
1098%                            TotalSup2Para(BoundarySup)*3-1;...
1099%                            TotalSup2Para(BoundarySup)*3];
1100%      BoundarySupPtrSub =BoundarySupPtrSub(:);
1101%      SubSupPtrSub = [ TotalSup2Para(SubSup)*3-2;...
1102%                       TotalSup2Para(SubSup)*3-1;...
1103%                       TotalSup2Para(SubSup)*3];
1104%      SubSupPtrSub =SubSupPtrSub(:);
1105   
1106      % clearn RayAllM PosiM CoPM1 HoriStickM_i VertStickM_i
1107%      NewRayAllM = RayAllM(:,TotalSupPtr);
1108      NewRayAllM = RayAllM(:,SubSupPtr);
1109      tar = sum(NewRayAllM ~= 0,2) ==3;
1110      NewRayAllM = NewRayAllM(tar,:);
1111
1112%      NewPosiMPPCP = PosiMPPCP(:,TotalSupPtr);
1113      NewPosiMPPCP = PosiMPPCP(:,SubSupPtr);
1114      tar = sum(NewPosiMPPCP ~= 0,2) ==3;
1115      NewPosiMPPCP = NewPosiMPPCP(tar,:);
1116      NewVarM = VarM2(tar);
1117
1118%      NewCoPM = CoPM1(:,TotalSupPtr) - CoPM2(:,TotalSupPtr);
1119      NewCoPM = CoPM1(:,SubSupPtr) - CoPM2(:,SubSupPtr);
1120      NewCoPMBound = CoPM1(:,BoundarySupPtr) - CoPM2(:,BoundarySupPtr);
1121%      tar = sum( NewCoPM ~= 0,2)  ==6;
1122      tar = sum( NewCoPM ~= 0,2) + sum( NewCoPMBound ~= 0,2) ==6;
1123      NewCoPM = NewCoPM(tar,:);
1124      NewCoPMBound = NewCoPMBound*PlanePara(BoundarySupPtr); % column vertor
1125      NewCoPMBound = NewCoPMBound(tar);
1126      NewCoPEstDepth = CoPEstDepth(tar);   
1127
1128%      NewHoriStickM = HoriStickM_i(:,TotalSupPtr)-HoriStickM_j(:,TotalSupPtr);
1129      NewHoriStickM = HoriStickM_i(:,SubSupPtr)-HoriStickM_j(:,SubSupPtr);
1130      NewHoriStickMBound = HoriStickM_i(:,BoundarySupPtr)-HoriStickM_j(:,BoundarySupPtr);
1131%      tar = sum(NewHoriStickM ~= 0,2) ==6;
1132      tar = sum(NewHoriStickM ~= 0,2) + sum( NewHoriStickMBound ~= 0,2)==6;
1133      NewHoriStickM = NewHoriStickM(tar,:);
1134      NewHoriStickMBound = NewHoriStickMBound*PlanePara(BoundarySupPtr); % column vertor
1135      NewHoriStickMBound = NewHoriStickMBound(tar);
1136      NewEstDepHoriStick = EstDepHoriStick(tar);
1137         NewWeightHoriNeighborStitch = WeightHoriNeighborStitch(tar);
1138
1139%      NewVertStickM = VertStickM_i(:,TotalSupPtr)-VertStickM_j(:,TotalSupPtr);
1140      NewVertStickM = VertStickM_i(:,SubSupPtr)-VertStickM_j(:,SubSupPtr);
1141      NewVertStickMBound = VertStickM_i(:, BoundarySupPtr)-VertStickM_j(:,BoundarySupPtr);
1142%      tar = sum(NewVertStickM ~= 0,2) ==6;
1143      tar = sum(NewVertStickM ~= 0,2) + sum(NewVertStickMBound ~= 0,2)==6;
1144      NewVertStickM = NewVertStickM(tar,:);
1145      NewVertStickMBound = NewVertStickMBound*PlanePara(BoundarySupPtr); % column vertor
1146      NewVertStickMBound = NewVertStickMBound(tar);
1147      NewEstDepVertStick = EstDepVertStick(tar);
1148         NewWeightVertNeighborStitch = WeightVertNeighborStitch(tar);
1149
1150%    try reduce the vertical constrain
1151     NonVertPtr = setdiff( 1:(3*NuSubSupSize), SubSup2Para( NuSup( (intersect( indexVertical,SubSupPtr ) +1)/3))*3-1 );
1152     YNoComingBack = YPointer(Sup2Para(SubSup));
1153     YNoComingBack(SubSup2Para( NuSup( (intersect( indexVertical,SubSupPtr ) +1)/3))) = [];
1154     YCompMask = zeros(1,3*NuSubSupSize);
1155     YCompMask(3*(1:NuSubSupSize)-1) = 1;
1156     YCompMask = YCompMask(NonVertPtr);
1157     XCompMask = zeros(1,3*NuSubSupSize);
1158     XCompMask( SubSup2Para( NuSup( (intersect( indexGroundX,SubSupPtr ) +2)/3))*3-2 ) = 1;
1159     XCompMask = XCompMask(NonVertPtr);
1160     ZCompMask = zeros(1,3*NuSubSupSize);
1161     ZCompMask( SubSup2Para( NuSup( (intersect( indexGroundZ,SubSupPtr ) )/3))*3 ) = 1;
1162     ZCompMask = ZCompMask(NonVertPtr);
1163     GroundMask = intersect( find(groundPara), ( SubSupPtr( 3:3:( size(SubSupPtr,1)) )/3 ) ) ;
1164
1165      Para = sdpvar( size(NonVertPtr',1),1);
1166%      F = set(Para(3*(1:NuSubSupSize)-1).*YPointer(Sup2Para(SubSup))<=0)...
1167      F = set(Para( logical(YCompMask)).*YNoComingBack<=0)...
1168          +set(NewRayAllM(:, NonVertPtr)*Para <=1/ClosestDist)...
1169          +set(NewRayAllM(:, NonVertPtr)*Para >=1/FarestDist);
1170%          +set( Para( SubSup2Para( NuSup( (intersect( indexVertical,SubSupPtr ) +1)/3))*3-1) == 0);
1171%          +set(Para(BoundarySupPtrSub) == PlanePara(BoundarySupPtr) )...
1172          % Special constrain for decomp fix the solved neighbor plane parameter       
1173
1174          %      +set(Para(indexGroundX) == 0)...
1175          %      +set(Para(indexGroundZ) == 0);
1176
1177%       if FractionalDepthError
1178% %         size(PosiMPPCP)
1179
1180    fprintf([ '     ' num2str( toc(inferenceTime) ) '\n     : In 2nd level Optimization, using ' opt.solver '.  ' ...
1181                     '(' num2str(k+1) '/' num2str((TotalRectX-1)+1) ',' num2str(l+1) '/' num2str((TotalRectY-1)+1) ')']);
1182       
1183    sol=solvesdp(F,norm( ( NewPosiMPPCP(:, NonVertPtr)*Para-ones(size(NewPosiMPPCP,1),1))./exp(abs(NewVarM)/BandWith),1)...
1184                 +Center*norm(( NewCoPM(:, NonVertPtr)*Para + NewCoPMBound).*NewCoPEstDepth, 1)...
1185                 +norm(( NewHoriStickM(:, NonVertPtr)*Para + NewHoriStickMBound).*...
1186                         NewEstDepHoriStick.*NewWeightHoriNeighborStitch,1)...
1187                 +norm(( NewVertStickM(:, NonVertPtr)*Para + NewVertStickMBound).*...
1188                         NewEstDepVertStick.*NewWeightVertNeighborStitch,1)...
1189                 +10*norm( ( Para( logical(XCompMask) ))./...
1190                            normPara( GroundMask )', 1)...
1191                 +10*norm( ( Para( logical(ZCompMask)))./...
1192                            normPara( GroundMask )', 1) ...
1193                 , opt);
1194%                 +10*norm( ( Para( XCompMask ))./...
1195%                            normPara( intersect( find(groundPara), ( SubSupPtr( 3:3:( size(SubSupPtr,1)) )/3 ) ) )', 1)...
1196%                 +10*norm( ( Para( ZCompMask))./...
1197%                            normPara( intersect( find(groundPara), ( SubSupPtr( 3:3:( size(SubSupPtr,1)) )/3 ) ) )', 1) ...
1198
1199             
1200             
1201             %             sqrt( max(CoPM1*ParaPPCP(:), 1/FarestDist).*max(CoPM2*ParaPPCP(:), 1/FarestDist)), 1)...
1202%             +norm(((VertStickM_i-VertStickM_j)*Para)./...
1203%             sqrt( max(VertStickM_i*ParaPPCP(:), 1/FarestDist).*max(VertStickM_j*ParaPPCP(:), 1/FarestDist)),1)...
1204%   pause;
1205    %         +1000*norm( Para(indexVertical)./ normPara( find(verticalPara) )',1) ...
1206%            +norm(((HoriStickM_i-HoriStickM_j)*Para)./sqrt((VertStickM_i*ParaPPCP(:)).*(VertStickM_j*ParaPPCP(:))),1)...
1207% %       else % not used anymore
1208%          solvesdp(F,norm( RayMd*Para - DepthInverseM,1)...
1209%                  +Center*norm((CoPM1 - CoPM2)*Para,1)...
1210%                  +norm((VertStickM_i-VertStickM_j)*Para,1)...
1211%                  +norm((HoriStickM_i-HoriStickM_j)*Para,1)...
1212%                  +1000*norm( Para(indexVertical)./ normPara( find(verticalPara) )',1) ...
1213%                  +1000*norm( Para(indexGroundX)./ normPara( find(groundPara) )', 1)  ...
1214%                  +1000*norm( Para(indexGroundZ)./ normPara( find(groundPara) )', 1) ...
1215%                  , opt);
1216% %         +0.01*norm( RayMtilt*ParaPPCP,1)...       
1217% %         +norm(spdiags(WeiCoP,0,size(CoPM1,1),size(CoPM1,1))*(CoPM1 - CoPM2)*ParaPPCP,1)...
1218% %         +norm( (CoPM1 - CoPM2)*ParaPPCP,1 )...
1219%      end
1220
1221     Para = double(Para);
1222     %sum(isnan(Para))
1223     yalmip('clear');
1224     tempPara = zeros(3*NuSubSupSize,1);
1225     tempPara(NonVertPtr) = Para;
1226%     PlanePara(TotalSupPtr) = Para;
1227     PlanePara(SubSupPtr) = tempPara;
1228    end
1229  end
1230
1231  %pause
1232  PlanePara = reshape(PlanePara,3,[]);
1233  FitDepth = FarestDist*ones(1,Default.VertYNuDepth*Default.HoriXNuDepth);
1234  FitDepth(~maskSkyEroded) = (1./sum(PlanePara(:,Sup2Para(SupEpand(~maskSkyEroded))).*Ray(:,~maskSkyEroded),1))';
1235  FitDepth = reshape(FitDepth,Default.VertYNuDepth,[]);
1236 
1237  % ==========Storage ==============
1238  if Default.Flag.AfterInferenceStorage
1239     save([ Default.ScratchFolder '/' strrep( Default.filename{1},'.jpg','') '_AInf.mat' ], 'FitDepth', 'depthMap', ...
1240            'Sup', 'SupOri', 'MedSup', 'RayOri','Ray','SupNeighborTable','maskSky','maskG','MultiScaleSupTable');
1241  end
1242  % ===============================
1243  %sum(isnan(FitDepth(:)))
1244  [Position3DFited] = im_cr2w_cr(FitDepth,permute(Ray,[2 3 1]));
1245  Position3DFited(3,:) = -Position3DFited(3,:);
1246  Position3DFited = permute(Position3DFited,[2 3 1]);
1247  RR =permute(Ray,[2 3 1]);
1248  temp = RR(:,:,1:2)./repmat(RR(:,:,3),[1 1 2]);
1249  PositionTex = permute(temp./repmat(cat(3,Default.a_default,Default.b_default),...
1250                [Default.VertYNuDepth Default.HoriXNuDepth 1])...
1251                +repmat(cat(3,Default.Ox_default,Default.Oy_default),...
1252                [Default.VertYNuDepth Default.HoriXNuDepth 1]),[3 1 2]);
1253  PositionTex = permute(PositionTex,[2 3 1]);
1254
1255
1256  fprintf(['     ' num2str( toc(inferenceTime) ) '\n     : Writing WRL.']);
1257  WrlFacestHroiReduce(Position3DFited, PositionTex, SupOri, Default.filename{1}, Default.filename{1}, ...
1258      Default.OutPutFolder, GridFlag, 0);
1259 
1260  system(['gzip -9 -c ' Default.OutPutFolder Default.filename{1} '.wrl > ' ...
1261         Default.OutPutFolder Default.filename{1} '.wrl.gz']);
1262  system(['cp '  Default.OutPutFolder Default.filename{1} '.wrl.gz ' ...
1263          Default.OutPutFolder Default.filename{1} '.wrl']);
1264  delete([Default.OutPutFolder Default.filename{1} '.wrl.gz']);
1265
1266    % Trial for segmentaion of the plane parameters
1267    if PlaneParaSegmentationFlag
1268       TempSup = Sup;%imresize(Sup,[330 610]);
1269       TempSup(TempSup==0) = max(TempSup(:))+1;
1270       TempSup2Para = Sup2Para;
1271       TempSup2Para(end+1) = size(PlanePara,2);
1272       TempPlanePara = PlanePara;
1273       TempPlanePara(:,end+1) = 0;
1274       PlaneParaPics = TempPlanePara(:,TempSup2Para( TempSup));
1275       PlaneParaPics = PlaneParaPics./repmat( 2*max(abs(PlaneParaPics),[],2), 1, size(PlaneParaPics,2));
1276       PlaneParaPics2 = permute( reshape( PlaneParaPics, 3, 55, []), [2 3 1])+0.5;
1277       MergedPlaneParaPics = segmentImgOpt( Default.sigm*scale, Default.k*scale, Default.minp*10*scale, PlaneParaPics, '', 0) + 1;
1278       if Default.Flag.DisplayFlag
1279          figure(400);
1280          imaegsc( PlaneParaPics);
1281          imaegsc( MergedPlaneParaPics);
1282       end
1283       save('/afs/cs/group/reconstruction3d/scratch/testE/PlaneParaPics2.mat','PlaneParaPics2','MergedPlaneParaPics');
1284    end
1285% save depth Map ++++++++++++++++
1286%       depthMap = FitDepth;
1287%       system(['mkdir ' Default.ScratchDataFolder '/VNonSupport_' DepthFolder '/']);
1288%       save([Default.ScratchDataFolder '/VNonSupport_' DepthFolder '/' depthfile '.mat'],'depthMap');
1289% =============================
1290% return;
Note: See TracBrowser for help on using the repository browser.