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

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

Added original make3d
File size: 23.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  ReportPlaneParaMRF(step, DepthFolder,...
40          Sup,MedSup,depthMap,VarMap,RayOri, Ray,MedRay,maskSky,maskG,Algo,k, CornerList, OccluList,...
41          MultiScaleSupTable, StraightLineTable, HBrokeBook, VBrokeBook,previoslyStored,...
42          baseline);
43% This function runs the RMF over the plane parameter of each superpixels
44
45global GeneralDataFolder ScratchDataFolder LocalFolder ClusterExecutionDirectory...
46    ImgFolder VertYNuPatch VertYNuDepth HoriXNuPatch HoriXNuDepth a_default b_default Ox_default Oy_default...
47    Horizon_default filename batchSize NuRow_default SegVertYSize SegHoriXSize WeiBatchSize PopUpVertY PopUpHoriX taskName;
48
49if nargin <20
50   baseline = 0;
51end
52
53% initialize parameters
54displayFlag = false;
55RenderVrmlFlag = true;
56StickHori = 5%0.1; % sticking power in horizontal direction
57StickVert = 5;     % sticking power in vertical direction
58Center = 5; % Co-Planar weight at the Center of each superpixel
59HoriConf = 1; % set the confidant of the learned depth at the middle in Horizontal direction of the image
60VertConf = 1; % set the confidant of the learned depth at the top of the image
61mapVert = linspace(VertConf,1,VertYNuDepth); % modeling the gravity prior
62mapHori = [linspace(HoriConf,1,round(HoriXNuDepth/2)) fliplr(linspace(HoriConf,1,HoriXNuDepth-round(HoriXNuDepth/2)))]; % modeling the gravity prior
63% ========set the range of depth that our model in
64ClosestDist = 1;
65% set the FarestDist to very 5 times to median depth
66FarestDist = 5*median(depthMap(:))
67% ================================================
68ceiling = 0*VertYNuDepth; % set the position of the ceiling, related to No plane coming back constrain
69Name{1} = 'FraWOPri';
70Name{2} = 'FraCoP';
71if isempty(MultiScaleSupTable)
72   Name{3} = 'FraStickCoP';
73else
74   Name{3} = 'FraStickCoPSTasCoP';
75end
76if ~isempty(MultiScaleSupTable)
77   MultiScaleFlag = true;
78   WeiV = ones(1,size(MultiScaleSupTable,2)-1);
79else
80   MultiScaleFlag = false;
81   WeiV = 1;
82end
83WeiV(1,1:2:end) = 1/3; % emphasize the middle scale three times smaller than large scale
84WeiV =WeiV./sum(WeiV);% normalize if pair of superpixels have same index in all the scale, their weight will be 10
85gravity =true; % if true, apply the HoriConf and VertConf linear scale weight
86CoPST = true; % if true, apply the Straight line prior as the Co-Planar constrain
87ConerImprove = false;
88FractionalDepthError = true;
89% get rid of the error point and sky point in the depthMap
90% set every depth bigger than FarestDistmeter to FarestDistmeters
91%CleanedDepthMap = (depthMapif ~previoslyStored >80).*medfilt2(depthMap,[4 4])+(depthMap<=80).*depthMap;
92CleanedDepthMap = depthMap;
93%CleanedDepthMap(depthMap>FarestDist) = FarestDist; % don't clean the point >80 sometimes it occlusion
94CleanedDepthMap(depthMap>FarestDist) = NaN; % don't clean the point >80 sometimes it occlusion
95Posi3D = im_cr2w_cr(CleanedDepthMap,permute(Ray,[2 3 1]));
96if ~previoslyStored
97
98NewMap = [rand(max(Sup(:)),3); [0 0 0]];
99% Clean the Sup near sky
100%maskSky = imdilate(maskSky, strel('disk', 3) );
101%[Sup,MedSup]=CleanedSup(Sup,MedSup,maskSky);
102maskSky = Sup == 0;
103maskSkyEroded = imerode(maskSky, strel('disk', 4) );
104SupEpand = ExpandSup2Sky(Sup,maskSkyEroded);
105NuPatch = HoriXNuDepth*VertYNuDepth-sum(maskSky(:));
106
107NuSup = setdiff(unique(Sup)',0);
108NuSup = sort(NuSup);
109NuSupSize = size(NuSup,2);
110
111% Sup index and planeParameter index inverse map
112Sup2Para = sparse(1,max(Sup(:)));
113Sup2Para(NuSup) = 1:NuSupSize;
114
115% constructinf the Straight line prior matrix Will be add in the CoPlane matrix
116NuLine = size(StraightLineTable,2);
117CoPSTList = [];
118for i = 1:NuLine
119    L = StraightLineTable{i};
120    X = L(1:(end-1))';
121    Y = L(2:end)';
122    if isempty(X)
123       continue;
124    end
125    for j = 1:size(X,1)
126        if X(j)~=Y(j)
127           CoPSTList = [CoPSTList; X(j) Y(j)];
128        end
129    end
130end
131end
132
133% Generate the Matrix for MRF
134tic
135PosiM = sparse(0,0);
136VarM = sparse(0,0);
137RayMd = sparse(0,0);
138RayAllOriM = sparse(0,0);
139RayAllM = sparse(0,0);
140RayMtilt = sparse(0,0);
141RayMCent = sparse(0,0);
142DepthInverseMCent = [];
143DepthInverseM = [];
144YPointer = [];
145beta = [];
146EmptyIndex = [];
147for i = NuSup
148%    mask = Sup ==i;
149    mask = SupEpand ==i; % include the Ray that will be use to expand the NonSky
150    RayAllOriM = blkdiag( RayAllOriM, RayOri(:,mask)');
151    RayAllM = blkdiag( RayAllM, Ray(:,mask)');
152    mask = Sup ==i; % Not include the Ray that will be use to expand the NonSky   
153    [yt x] = find(mask);
154    CenterX = round(median(x));
155    CenterY = round(median(yt));
156    YPointer = [YPointer; CenterY >= ceiling];
157    mask(isnan(CleanedDepthMap)) = false;
158    SupNuPatch(i) = sum(mask(:));
159%    if sum(mask(:)) < 5
160%       EmptyIndex = [EmptyIndex; i];
161%       mask(:) = false;
162%    end
163    % find center point
164    [yt x] = find(mask);
165    CenterX = round(median(x));
166    CenterY = round(median(yt));
167 
168  if ~all(mask(:)==0)
169    if gravity
170      Pa2CenterRatio = median(CleanedDepthMap(mask))./CleanedDepthMap(mask);
171      if sum(mask(:)) > 0
172         RayMtilt = blkdiag(RayMtilt, ...
173             ( Pa2CenterRatio(:,[1 1 1]).*repmat(RayOri(:,CenterY,CenterX)',[ SupNuPatch(i) 1])- RayOri(:,mask)'));
174      else
175         RayMtilt = blkdiag(RayMtilt, RayOri(:,mask)');
176      end
177      RayMCent = blkdiag(RayMCent, RayOri(:,CenterY,CenterX)'*SupNuPatch(i)*mapVert(CenterY)*mapHori(CenterX));
178      PosiM = blkdiag(PosiM,Posi3D(:,mask)'.*repmat( mapVert(yt)',[1 3]).*repmat( mapHori(x)',[1 3]));
179      VarM = [VarM; VarMap(mask)];
180      RayMd = blkdiag(RayMd,RayOri(:,mask)'.*repmat( mapVert(yt)',[1 3]).*repmat( mapHori(x)',[1 3]));
181      DepthInverseMCent = [DepthInverseMCent; 1./median(CleanedDepthMap(mask)).*SupNuPatch(i)* mapVert(CenterY)'*mapHori(CenterX)];
182      DepthInverseM = [DepthInverseM; 1./CleanedDepthMap(mask).* mapVert(yt)'.*mapHori(x)'];
183    else
184      Pa2CenterRatio = CleanedDepthMap(CenterY,CenterX)./CleanedDepthMap(mask);
185      if sum(mask(:)) > 0
186         RayMtilt = blkdiag(RayMtilt, ...
187             ( Pa2CenterRatio(:,[1 1 1]).*repmat(RayOri(:,CenterY,CenterX)',[ SupNuPatch(i) 1])- RayOri(:,mask)'));
188      alse
189         RayMtilt = blkdiag(RayMtilt, RayOri(:,mask)');
190      end
191      RayMCent = blkdiag(RayMCent, RayOri(:,CenterY,CenterX)'*SupNuPatch(i));
192      PosiM = blkdiag(PosiM,Posi3D(:,mask)');
193      VarM = [VarM; VarMap(mask)];
194      RayMd = blkdiag(RayMd,RayOri(:,mask)');
195      DepthInverseMCent = [DepthInverseMCent; 1./median(CleanedDepthMap(mask)).*SupNuPatch(i)];
196      DepthInverseM = [DepthInverseM; 1./CleanedDepthMap(mask)];
197    end
198  else
199     RayMtilt = blkdiag(RayMtilt, RayOri(:,mask)');
200     RayMCent = blkdiag(RayMCent, RayOri(:,mask)');
201     PosiM = blkdiag(PosiM, Posi3D(:,mask)');
202     VarM = [VarM; VarMap(mask)];
203     RayMd = blkdiag(RayMd, RayOri(:,mask)');
204  end
205end
206YPointer(YPointer==0) = -1;
207
208% buliding CoPlane Matrix=========================================================================
209BounaryPHori = conv2(Sup,[1 -1],'same') ~=0;
210BounaryPHori(:,end) = 0;
211BounaryPVert = conv2(Sup,[1; -1],'same') ~=0;
212BounaryPVert(end,:) = 0;
213ClosestNList = [ Sup(find(BounaryPHori==1)) Sup(find(BounaryPHori==1)+VertYNuDepth);...
214                 Sup(find(BounaryPVert==1)) Sup(find(BounaryPVert==1)+1)];
215ClosestNList = sort(ClosestNList,2);
216ClosestNList = unique(ClosestNList,'rows');
217ClosestNList(ClosestNList(:,1) == 0,:) = [];
218BoundaryAll = BounaryPHori + BounaryPHori(:,[1 1:(end-1)])...
219             +BounaryPVert + BounaryPVert([1 1:(end-1)],:);
220BoundaryAll([1 end],:) = 1;
221BoundaryAll(:,[1 end]) = 1;
222NuSTList = 0;
223if CoPST
224   ClosestNList = [ClosestNList; CoPSTList];
225   NuSTList = size(CoPSTList,1)
226end
227NuNei = size(ClosestNList,1);
228CoPM1 = sparse(0,3*NuSupSize);
229CoPM2 = sparse(0,3*NuSupSize);
230CoPEstDepth = sparse(0,0);
231CoPNorM = [];
232WeiCoP = [];
233for i = 1: NuNei
234%  if ~CornerList(i)
235    mask = Sup == ClosestNList(i,1);
236    SizeMaskAll = sum(mask(:));
237    [y x] = find(mask);
238    CenterX = round(median(x));
239    CenterY = round(median(y));
240    y = find(mask(:,CenterX));
241    if ~isempty(y)
242       CenterY = round(median(y));
243    end
244%    CenterX = round(median(x));
245%    CenterY = round(median(yt));
246   
247    temp1 = sparse(1, 3*NuSupSize);
248    temp2 = sparse(1, 3*NuSupSize);
249    temp1(:,(Sup2Para(ClosestNList(i,1))*3-2): Sup2Para(ClosestNList(i,1))*3) = Ray(:,CenterY,CenterX)';
250    temp2(:,(Sup2Para(ClosestNList(i,2))*3-2): Sup2Para(ClosestNList(i,2))*3) = Ray(:,CenterY,CenterX)';
251    NuNei-NuSTList
252    i
253    if i < NuNei-NuSTList % only immediate connecting superpixels are weighted using MultiScaleSup
254%       wei = WeiV*10;%*(MultiScaleSupTable(Sup2Para(ClosestNList(i,1)),2:end) == MultiScaleSupTable(Sup2Para(ClosestNList(i,2)),2:end))'; 
255       if MultiScaleFlag
256          wei = WeiV*(MultiScaleSupTable(Sup2Para(ClosestNList(i,1)),2:end) == MultiScaleSupTable(Sup2Para(ClosestNList(i,2)),2:end))'
257       else
258          wei = WeiV
259       end
260    else
261       wei = sum(WeiV)
262    end
263    CoPM1 = [CoPM1; temp1*wei];
264    CoPM2 = [CoPM2; temp2*wei];
265    tempWeiCoP = [SizeMaskAll];
266    CoPEstDepth = [CoPEstDepth; max(median(CleanedDepthMap(mask)),ClosestDist)];
267   
268    mask = Sup == ClosestNList(i,2);
269    SizeMaskAll = sum(mask(:));
270    [y x] = find(mask);
271    CenterX = round(median(x));
272    CenterY = round(median(y));
273    y = find(mask(:,CenterX));
274    if ~isempty(y)
275       CenterY = round(median(y));
276    end
277%    [yt x] = find(mask);
278%    CenterX = round(median(x));
279%    CenterY = round(median(yt));
280
281    temp1 = sparse(1, 3*NuSupSize);
282    temp2 = sparse(1, 3*NuSupSize);
283    temp1(:,(Sup2Para(ClosestNList(i,1))*3-2): Sup2Para(ClosestNList(i,1))*3) = Ray(:,CenterY,CenterX)';
284    temp2(:,(Sup2Para(ClosestNList(i,2))*3-2): Sup2Para(ClosestNList(i,2))*3) = Ray(:,CenterY,CenterX)';
285    CoPM1 = [CoPM1; temp1*wei];
286    CoPM2 = [CoPM2; temp2*wei];
287    tempWeiCoP = [tempWeiCoP; SizeMaskAll];
288    WeiCoP = [WeiCoP; tempWeiCoP];
289    CoPEstDepth = [CoPEstDepth; max(median(CleanedDepthMap(mask)),ClosestDist)];
290%  end
291end%=========================================================================================================
292
293% find the boundary point that might need to be stick ot each other==========================================
294HoriStickM_i = sparse(0,3*NuSupSize);
295HoriStickM_j = sparse(0,3*NuSupSize);
296HoriStickPointInd = [];
297EstDepHoriStick = [];
298for i = find(BounaryPHori==1)'
299    j = i+VertYNuDepth;
300    if Sup(i) == 0 || Sup(j) == 0
301       continue;
302    end
303%  if ~OccluList(sum( ClosestNList == repmat(sort([Sup(i) Sup(j)]), [NuNei  1]),2) == 2)
304%  size(OccluList)
305%  if ~any(sum( OccluList == repmat(sort([Sup(i) Sup(j)]), [size(OccluList,1)  1]),2) == 2)
306    Target(1) = Sup2Para(Sup(i));
307    Target(2) = Sup2Para(Sup(j));
308    rayBoundary(:,1) =  Ray(:,i);
309    rayBoundary(:,2) =  Ray(:,i);
310%    betaTemp = StickHori;%*(DistStickLengthNormWei.^2)*beta(Target(I));
311%    betaTemp = StickHori*WeiV;%*(MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end))';%*(DistStickLengthNormWei.^2)*beta(Target(I));
312    if MultiScaleFlag
313       betaTemp = StickHori*WeiV*(MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end))';%*(DistStickLengthNormWei.^2)*beta(Target(I));
314    else
315       betaTemp = StickHori;
316    end
317    temp = sparse(3,NuSupSize);
318    temp(:,Target(1)) = rayBoundary(:,1);
319    HoriStickM_i = [HoriStickM_i; betaTemp*temp(:)'];
320    temp = sparse(3,NuSupSize);
321    temp(:,Target(2)) = rayBoundary(:,2);
322    HoriStickM_j = [HoriStickM_j; betaTemp*temp(:)'];
323    EstDepHoriStick = [EstDepHoriStick; sqrt(max(CleanedDepthMap(i),ClosestDist)*max(CleanedDepthMap(j),ClosestDist))];
324    HoriStickPointInd = [HoriStickPointInd i ];
325%  else
326%    disp('Occlu');
327%  end
328end
329VertStickM_i = sparse(0,3*NuSupSize);
330VertStickM_j = sparse(0,3*NuSupSize);
331VertStickPointInd = [];
332EstDepVertStick = [];
333for i = find(BounaryPVert==1)'
334    j = i+1;
335    if Sup(i) == 0 || Sup(j) == 0
336       continue;
337    end
338%  if ~OccluList(sum( ClosestNList == repmat(sort([Sup(i) Sup(j)]), [NuNei  1]),2) == 2)
339%  if ~any(sum( OccluList == repmat(sort([Sup(i) Sup(j)]), [size(OccluList,1)  1]),2) == 2)
340    Target(1) = Sup2Para(Sup(i));
341    Target(2) = Sup2Para(Sup(j));
342    rayBoundary(:,1) =  Ray(:,i);
343    rayBoundary(:,2) =  Ray(:,i);
344%    betaTemp = StickVert;%*(DistStickLengthNormWei.^2)*beta(Target(I));
345    if MultiScaleFlag
346       betaTemp = StickVert*WeiV*(MultiScaleSupTable(Sup2Para(Sup(i)),2:end) == MultiScaleSupTable(Sup2Para(Sup(j)),2:end))';
347    else
348       betaTemp = StickVert;
349    end
350    temp = sparse(3,NuSupSize);
351    temp(:,Target(1)) = rayBoundary(:,1);
352    VertStickM_i = [VertStickM_i; betaTemp*temp(:)'];
353    temp = sparse(3,NuSupSize);
354    temp(:,Target(2)) = rayBoundary(:,2);
355    VertStickM_j = [VertStickM_j; betaTemp*temp(:)'];
356    EstDepVertStick = [EstDepVertStick; sqrt(max(CleanedDepthMap(i),ClosestDist)*max(CleanedDepthMap(j),ClosestDist))];
357    VertStickPointInd = [VertStickPointInd i ];
358%  else
359%    disp('Occlu');
360%  end
361end
362% ======================================Finish building up matrix=====================hard work======================
363% ================================================================================================================
364depthfile = strrep(filename{k},'img','depth_learned'); %
365
366for i = step
367    opt = sdpsettings('solver','sedumi');
368    ParaPPCP = sdpvar(3*NuSupSize,1);
369    F = set(ParaPPCP(3*(1:NuSupSize)-1).*YPointer<=0) ...
370        +set(RayAllM*ParaPPCP >=1/FarestDist)...
371        +set(RayAllM*ParaPPCP <=1/ClosestDist)...
372        +set(RayAllOriM*ParaPPCP >=1/FarestDist)...
373        +set(RayAllOriM*ParaPPCP <=1/ClosestDist);
374% First fit the plane to find the estimated plane parameters
375% If PosiM contain NaN data the corresponding Plane Parameter will also be NaN
376    if i == 1 % W/O  priors
377    solvesdp(F,norm( PosiM*ParaPPCP-ones(size(PosiM,1),1),1)./(VarMap)...
378             , opt);
379    elseif i ==2 % Coner
380    solvesdp(F,norm( PosiM*ParaPPCP-ones(size(PosiM,1),1),1)./(VarMap)...
381             +Center*norm(((CoPM1 - CoPM2)*ParaPPCP).*CoPEstDepth, 1)...
382             , opt);
383    else % all the priors
384    solvesdp(F,norm( PosiM*ParaPPCP-ones(size(PosiM,1),1),1)./(VarMap)...
385             +Center*norm(((CoPM1 - CoPM2)*ParaPPCP).*CoPEstDepth, 1)...
386             +norm(((HoriStickM_i-HoriStickM_j)*ParaPPCP).*EstDepHoriStick,1)...
387             +norm(((VertStickM_i-VertStickM_j)*ParaPPCP).*EstDepVertStick,1)...
388             , opt);
389    end
390
391    ParaPPCP = double(ParaPPCP);
392    SepPointMeasureHori = 1./(HoriStickM_i*ParaPPCP)-1./(HoriStickM_j*ParaPPCP);
393    SepPointMeasureVert = 1./(VertStickM_i*ParaPPCP)-1./(VertStickM_j*ParaPPCP);
394    ParaPPCP = reshape(ParaPPCP,3,[]);
395    %any(any(isnan(ParaPPCP)))
396    % porject the ray on planes to generate the ProjDepth
397    FitDepthPPCP = FarestDist*ones(1,VertYNuDepth*HoriXNuDepth);
398    FitDepthPPCP(~maskSkyEroded) = (1./sum(ParaPPCP(:,Sup2Para(SupEpand(~maskSkyEroded))).*Ray(:,~maskSkyEroded),1))';
399    FitDepthPPCP = reshape(FitDepthPPCP,VertYNuDepth,[]);
400    % storage for comparison purpose ======================
401    depthMap = FarestDist*ones(1,VertYNuDepth*HoriXNuDepth);
402    depthMap(~maskSkyEroded) = (1./sum(ParaPPCP(:,Sup2Para(SupEpand(~maskSkyEroded))).*RayOri(:,~maskSkyEroded),1))';
403    depthMap = reshape(depthMap,VertYNuDepth,[]);
404%    if baseline == 0
405       system(['mkdir ' ScratchDataFolder '/' Name{i} '_' DepthFolder '/']);
406       save([ScratchDataFolder '/' Name{i} '_' DepthFolder '/' depthfile '.mat'],'depthMap');
407%    elseif baseline == 1
408%       system(['mkdir ' ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline/']);
409%       save([ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline/' depthfile '.mat'],'depthMap');
410%    else
411%       system(['mkdir ' ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline2/']);
412%       save([ScratchDataFolder '/' Name{i} '_' DepthFolder '_baseline2/' depthfile '.mat'],'depthMap');
413%    end
414 
415    % =====================================================
416    [Position3DFitedPPCP] = im_cr2w_cr(FitDepthPPCP,permute(Ray,[2 3 1]));
417    if RenderVrmlFlag && i==3
418       [VrmlName] = vrml_test_faceset_goodSkyBoundary(  filename{k}, Position3DFitedPPCP, FitDepthPPCP, permute(Ray,[2 3 1]), [Name{i} DepthFolder], ...
419                    [], [], 0, maskSkyEroded, maskG, 1, 0, a_default, b_default, Ox_default, Oy_default);
420       system(['gzip -9 -c ' ScratchDataFolder '/vrml/' VrmlName ' > ' ScratchDataFolder '/vrml/' VrmlName '.gz']);
421       delete([ScratchDataFolder '/vrml/' VrmlName]);
422    end
423%    save([ScratchDataFolder '/data/PreDepth/FitDepthPPCP' num2str(k) '.mat'],'FitDepthPPCP','SepPointMeasureHori',...
424%      'SepPointMeasureVert','VertStickPointInd','HoriStickPointInd');
425%save([ScratchDataFolder '/data/All.mat']);
426end
427
428%return;
429%==================Finished for one step MRF==========================================================================================================
430
431
432% ====================================following are 2nd and 3rd step MRF to give more visually pleasing result=======================================
433if ConerImprove
434[CoPM1 CoPM2 WeiCoP Sup Sup2ParaNew FixPara NuSup VertStickPointInd HoriStickPointInd] ...
435            = FixParaFreeCorner( Sup, Ray, HoriStickPointInd, VertStickPointInd,...
436              SepPointMeasureHori, SepPointMeasureVert, OccluList, WeiV);
437
438RayAllM = sparse(0,0);
439YPointer = [];
440for i = NuSup
441    mask = Sup ==i;
442    RayAllM = blkdiag( RayAllM, Ray(:,mask)');
443    [yt x] = find(mask);
444    CenterX = round(median(x));
445    CenterY = round(median(yt));
446    YPointer = [YPointer; CenterY >= ceiling];
447end
448
449NuSupSize = size(NuSup,2);
450opt = sdpsettings('solver','sedumi');
451ParaCorner = sdpvar(3*NuSupSize,1);
452F = set(ParaCorner(3*(1:NuSupSize)-1).*YPointer<=0) + set(RayAllM*ParaCorner >=1/FarestDist)...
453    +set(RayAllM*ParaCorner <=1/ClosestDist)...
454    +set(ParaCorner([(Sup2ParaNew(FixPara)*3-2); (Sup2ParaNew(FixPara)*3-1); (Sup2ParaNew(FixPara)*3)]) == ...
455     ParaPPCP([(Sup2Para(FixPara)*3-2); (Sup2Para(FixPara)*3-1); (Sup2Para(FixPara)*3)]));
456solvesdp(F,Center*norm((CoPM1 - CoPM2)*ParaCorner, 1)...
457          , opt);
458ParaCorner = double(ParaCorner);
459ParaCorner = reshape(ParaCorner,3,[]);
460FitDepthCorner = 80*ones(1,VertYNuDepth*HoriXNuDepth);
461FitDepthCorner(Sup~=0) = (1./sum(ParaCorner(:,Sup2ParaNew(Sup(Sup~=0))).*Ray(:,Sup~=0),1))';
462FitDepthCorner = reshape(FitDepthCorner,VertYNuDepth,[]);
463[Position3DFitedCorner] = im_cr2w_cr(FitDepthCorner,permute(Ray,[2 3 1]));
464[VrmlName] = vrml_test_faceset_goodSkyBoundary( filename{k}, Position3DFitedCorner, FitDepthCorner, permute(Ray,[2 3 1]), [Name 'Corner'], ...
465[], [], 0, maskSkyEroded, maskG, 1, 0, a_default, b_default, Ox_default, Oy_default);
466system(['gzip -9 -c ' ScratchDataFolder '/vrml/' VrmlName ' > ' ScratchDataFolder '/vrml/' VrmlName '.gz']);
467delete([ScratchDataFolder '/vrml/' VrmlName]);
468end
469%return;
470% ============================================================
471;% generating new PosiMPPCP using the new position
472PosiMPPCP = sparse(0,0);
473for i = NuSup
474    mask = Sup == i;
475    PosiMPPCP = blkdiag(PosiMPPCP, Position3DFitedPPCP(:,mask)');
476end
477
478
479  groundThreshold = cos(15*pi/180);
480  verticalThreshold = cos(50*pi/180);
481  normPara = norms(ParaPPCP);
482  normalizedPara = ParaPPCP ./ repmat( normPara, [3 1]);
483  groundPara = abs(normalizedPara(2,:)) >= groundThreshold;
484  verticalPara = abs(normalizedPara(2,:)) <= verticalThreshold;
485
486  indexVertical = find( verticalPara)*3-1;
487  indexGroundX = find( groundPara)*3-2;
488  indexGroundZ = find( groundPara)*3;
489  opt = sdpsettings('solver','sedumi');
490  Para = sdpvar(3*NuSupSize,1);
491  F = set(Para(3*(1:NuSupSize)-1).*YPointer<=0) + set(RayAllM*Para >=1/FarestDist)...
492      +set(RayAllM*Para <=1/ClosestDist);
493if FractionalDepthError
494   size(PosiMPPCP)
495   solvesdp(F,norm( PosiMPPCP*Para-ones(size(PosiMPPCP,1),1),1)...
496             +Center*norm(((CoPM1 - CoPM2)*Para)./...
497             sqrt( max(CoPM1*ParaPPCP(:), 1/FarestDist).*max(CoPM2*ParaPPCP(:), 1/FarestDist)), 1)...
498             +norm(((VertStickM_i-VertStickM_j)*Para)./...
499             sqrt( max(VertStickM_i*ParaPPCP(:), 1/FarestDist).*max(VertStickM_j*ParaPPCP(:), 1/FarestDist)),1)...
500             +1000*norm( Para(indexVertical)./ normPara( find(verticalPara) )',1) ...
501             +100*norm( Para(indexGroundX)./ normPara( find(groundPara) )', 1)  ...
502             +100*norm( Para(indexGroundZ)./ normPara( find(groundPara) )', 1) ...
503             , opt);
504%            +norm(((HoriStickM_i-HoriStickM_j)*Para)./sqrt((VertStickM_i*ParaPPCP(:)).*(VertStickM_j*ParaPPCP(:))),1)...
505else
506   solvesdp(F,norm( RayMd*Para - DepthInverseM,1)...
507            +Center*norm((CoPM1 - CoPM2)*Para,1)...
508            +norm((VertStickM_i-VertStickM_j)*Para,1)...
509            +norm((HoriStickM_i-HoriStickM_j)*Para,1)...
510            +1000*norm( Para(indexVertical)./ normPara( find(verticalPara) )',1) ...
511            +1000*norm( Para(indexGroundX)./ normPara( find(groundPara) )', 1)  ...
512            +1000*norm( Para(indexGroundZ)./ normPara( find(groundPara) )', 1) ...
513            , opt);
514%         +0.01*norm( RayMtilt*ParaPPCP,1)...       
515%         +norm(spdiags(WeiCoP,0,size(CoPM1,1),size(CoPM1,1))*(CoPM1 - CoPM2)*ParaPPCP,1)...
516%         +norm( (CoPM1 - CoPM2)*ParaPPCP,1 )...
517end
518  Para = reshape(Para,3,[]);
519  Para = double(Para);
520
521  FitDepth = 80*ones(1,VertYNuDepth*HoriXNuDepth);
522  FitDepth(~maskSkyEroded) = (1./sum(Para(:,Sup2Para(SupEpand(~maskSkyEroded))).*Ray(:,~maskSkyEroded),1))';
523  FitDepth = reshape(FitDepth,VertYNuDepth,[]);
524  sum(isnan(FitDepth(:)))
525  [Position3DFited] = im_cr2w_cr(FitDepth,permute(Ray,[2 3 1]));
526  [VrmlName] = vrml_test_faceset_goodSkyBoundary(       filename{k}, Position3DFited, FitDepth, permute(Ray,[2 3 1]), 'VH_WOHSTickMerged', ...
527                                        [], [], 0, maskSkyEroded, maskG, 1, 0, a_default, b_default, Ox_default, Oy_default);
528  system(['gzip -9 -c ' ScratchDataFolder '/vrml/' VrmlName ' > ' ScratchDataFolder '/vrml/' VrmlName '.gz']);
529  delete([ScratchDataFolder '/vrml/' VrmlName]);
530
531%end
532return;
Note: See TracBrowser for help on using the repository browser.