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

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

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