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

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

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