Hi,

I posted a message yesterday on the subject, but I didn't get any response, so I will try again to explain my problem. I'm sorry if I'm irritating anyone by posting two messages on the same problem, but since I got no answers, I think it's ok.

I'm currently programming a GUI for nivelling fields. Right now, my program draw the actual field using surf() and the best-fit plan using regress() and mesh(). The user can then draw a line on the best-fit plan to divide the field into 2 sections, so that the program calculate two best-fit planes for each of of the field section. 

My problem is, when the user draw the line, the GUI traces 2 other planes instead of the initial best-fit plane, but these planes are independant at their junctions, because I use the function regress() twice with differents parts of the field. I would like to do a 2-parts regression in 3d where the 2 best-fit planes have the same intersection line. Is there a function I can use for that? I read about the Statistic Toolbox, but I didn't found a function wich does Multi-phase/segmented mutiple linear regression in 3D.

I have also tried downloading some regression programs, like the mancovan program, wich include the mStepwise.m function, but I didn't succes at making it works.

I would very apreciate if somebody responds me and give me an idea on how could I program a segmented multiple linear regression, or the name of a function you know that can handle that, or anything else. Thank you.

Best regards,

P.S.: I'm sorry if I made some English mistakes writing this message. English is my second language.

"Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <hur1uk$qto$1@fred.mathworks.com>...
> Hi,
> 
> I posted a message yesterday on the subject, but I didn't get any response, so I will try again to explain my problem. I'm sorry if I'm irritating anyone by posting two messages on the same problem, but since I got no answers, I think it's ok.
> 
> I'm currently programming a GUI for nivelling fields. Right now, my program draw the actual field using surf() and the best-fit plan using regress() and mesh(). The user can then draw a line on the best-fit plan to divide the field into 2 sections, so that the program calculate two best-fit planes for each of of the field section. 
> 
> My problem is, when the user draw the line, the GUI traces 2 other planes instead of the initial best-fit plane, but these planes are independant at their junctions, because I use the function regress() twice with differents parts of the field. I would like to do a 2-parts regression in 3d where the 2 best-fit planes have the same intersection line. Is there a function I can use for that? I read about the Statistic Toolbox, but I didn't found a function wich does Multi-phase/segmented mutiple linear regression in 3D.
> 
> I have also tried downloading some regression programs, like the mancovan program, wich include the mStepwise.m function, but I didn't succes at making it works.
> 
> I would very apreciate if somebody responds me and give me an idea on how could I program a segmented multiple linear regression, or the name of a function you know that can handle that, or anything else. Thank you.
> 
> Best regards,
> 
> P.S.: I'm sorry if I made some English mistakes writing this message. English is my second language.

You can do the regression:

z ~ a1*x + b1*y + c1  for { (x,y,z) data on set # 1, e.g., left side of the line }
z ~ a2*x + b2*y + c2  for { (x,y,z) data on set # 2, e.g., right side of the line }

Under the constraints

a1*x + b1*y + c1 = a2*x + b2*y + c2 for { (x,y) two (arbitrary) points on the line }.

Bruno

> 
> You can do the regression:
> 
> z ~ a1*x + b1*y + c1  for { (x,y,z) data on set # 1, e.g., left side of the line }
> z ~ a2*x + b2*y + c2  for { (x,y,z) data on set # 2, e.g., right side of the line }
> 
> Under the constraints
> 
> a1*x + b1*y + c1 = a2*x + b2*y + c2 for { (x,y) two (arbitrary) points on the line }.
> 
> Bruno

Thank you. 
I understand what you say, because defining the constraint will force the coefficients to be dependent, so the regression planes will pass through the same line...

But how do I set this constraint?

According to my Matlab experience, I think I would do a FOR or WHILE loop to calculate the a1,b1,c1 and a2,b2,c2 coefficients a lots of time, each time evaluating the 'constraint' and see if it is respected, then finish the loop.

Is there another way to set this constraint?

Thank you for your help!

"Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <hur3j4$gno$1@fred.mathworks.com>...
> > 
> > You can do the regression:
> > 
> > z ~ a1*x + b1*y + c1  for { (x,y,z) data on set # 1, e.g., left side of the line }
> > z ~ a2*x + b2*y + c2  for { (x,y,z) data on set # 2, e.g., right side of the line }
> > 
> > Under the constraints
> > 
> > a1*x + b1*y + c1 = a2*x + b2*y + c2 for { (x,y) two (arbitrary) points on the line }.
> > 
> > Bruno
> 
> Thank you. 
> I understand what you say, because defining the constraint will force the coefficients to be dependent, so the regression planes will pass through the same line...
> 
> But how do I set this constraint?
> 
> According to my Matlab experience, I think I would do a FOR or WHILE loop to calculate the a1,b1,c1 and a2,b2,c2 coefficients a lots of time, each time evaluating the 'constraint' and see if it is respected, then finish the loop.
> 
> Is there another way to set this constraint?

Yes, just write down the linear system with constraints. The code  should goes like this (I have not checked it)

% (x,y,x) are three input vectors of data
% Assuming the line is x=y
set1 = x <= y; % upper
set2 = ~set1; % lower

x = x(:);
y = y(:);
z = z(:);

M1 = [x(set1) y(set1) ones(sum(set1),1)];
M2 = [x(set2) y(set2) ones(sum(set2),1)];
M = blkdiag(M1,M2);
Aeq = [0 0 1 0 0 -1; % match z at (x,y)=(0,0)
       1 1 1 -1 -1 -1]; % match z at (x,y)=(1,1)
B = null(Aeq);
rhs = [z(set1); z(set2)];
sol = B*((M*B)\z);
a1 = sol(1)
b1 = sol(2)
c1 = sol(3)
a2 = sol(4)
b2 = sol(5)
c2 = sol(6)

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hur5l9$4g9$1@fred.mathworks.com>...

> 
> Yes, just write down the linear system with constraints. The code  should goes like this (I have not checked it)
> 
> % (x,y,x) are three input vectors of data
> % Assuming the line is x=y
> set1 = x <= y; % upper
> set2 = ~set1; % lower
> 
> x = x(:);
> y = y(:);
> z = z(:);
> 
> M1 = [x(set1) y(set1) ones(sum(set1),1)];
> M2 = [x(set2) y(set2) ones(sum(set2),1)];
> M = blkdiag(M1,M2);
> Aeq = [0 0 1 0 0 -1; % match z at (x,y)=(0,0)
>        1 1 1 -1 -1 -1]; % match z at (x,y)=(1,1)
> B = null(Aeq);
> rhs = [z(set1); z(set2)];
> sol = B*((M*B)\z);
> a1 = sol(1)
> b1 = sol(2)
> c1 = sol(3)
> a2 = sol(4)
> b2 = sol(5)
> c2 = sol(6)
> 
> % Bruno

Thanks

I understand what you're doing, but in the division line where you solves the problem, shouldn't the last term be rhs instead of z? 
i.e. sol = B*((M*B)\rhs); instead of sol = B*((M*B)\z); ?

Also, I've never utilised the function null. Wouldn't the code work without this function?
i.e. simply write   sol = M\rhs;  ?

"Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <hur99g$6m7$1@fred.mathworks.com>...

> 
> I understand what you're doing, but in the division line where you solves the problem, shouldn't the last term be rhs instead of z? 
> i.e. sol = B*((M*B)\rhs); instead of sol = B*((M*B)\z); ?

Yes.

> 
> Also, I've never utilised the function null. Wouldn't the code work without this function?
> i.e. simply write   sol = M\rhs;  ?

No, because the constraints are not generally satisfied, so the two surface will not match et the line. That's what you are asking, no?

The constraints are satisfied if and only if A*sol = 0; i.e.,  sol belongs to span<B>.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hura1c$qus$1@fred.mathworks.com>...
> "Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <hur99g$6m7$1@fred.mathworks.com>...
> 
> > 
> > I understand what you're doing, but in the division line where you solves the problem, shouldn't the last term be rhs instead of z? 
> > i.e. sol = B*((M*B)\rhs); instead of sol = B*((M*B)\z); ?
> 
> Yes.
> 
> > 
> > Also, I've never utilised the function null. Wouldn't the code work without this function?
> > i.e. simply write   sol = M\rhs;  ?
> 
> No, because the constraints are not generally satisfied, so the two surface will not match et the line. That's what you are asking, no?
> 
> The constraints are satisfied if and only if A*sol = 0; i.e.,  sol belongs to span<B>.
> 
> Bruno

I tried the code with the GPS data of the field I have. It worked fine with the function null, but the plane were not matching at the line if the line i just sol = M\rhs;
As you said :P!

Since the possibility of separating a field into smaller sections is the main function of my field nivelling program, it helped me a lot.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hura1c$qus$1@fred.mathworks.com>...
> "Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <hur99g$6m7$1@fred.mathworks.com>...
> 
> > 
> > I understand what you're doing, but in the division line where you solves the problem, shouldn't the last term be rhs instead of z? 
> > i.e. sol = B*((M*B)\rhs); instead of sol = B*((M*B)\z); ?
> 
> Yes.
> 
> > 
> > Also, I've never utilised the function null. Wouldn't the code work without this function?
> > i.e. simply write   sol = M\rhs;  ?
> 
> No, because the constraints are not generally satisfied, so the two surface will not match et the line. That's what you are asking, no?
> 
> The constraints are satisfied if and only if A*sol = 0; i.e.,  sol belongs to span<B>.
> 
> Bruno

Hey Bruno, I have a question about your code.

when you have writen:
    Aeq = [0 0 1 0 0 -1; % match z at (x,y)=(0,0)
           1 1 1 -1 -1 -1]; % match z at (x,y)=(1,1)
    B = null(Aeq);
To define Aeq do you have to know that the equation of the line was x=y?

I believe so, but I want to developp is a code that works for any line. I know 2 points of the line, so I can obtain a vector that passes threw the line. After that, how can I write the matrix Aeq? I propose the following code, but I would like to know if it's correct, because I'm new to the function null so I don't understand all the concepts involved:
%I have p1(x,y,z) and p2(x,y,z) are two points from the line
%Can I write:
Aeq = [p1(1) p1(2) p1(3) -p1(1) -p1(2) -p1(3); % match z at p1?
           p2(1) p2(2) p2(3) -p2(1) -p2(2) -p2(3)];  % match z at p2?
%???

"Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <hurdtl$egg$1@fred.mathworks.com>...

> %Can I write:
> Aeq = [p1(1) p1(2) p1(3) -p1(1) -p1(2) -p1(3); % match z at p1?
>            p2(1) p2(2) p2(3) -p2(1) -p2(2) -p2(3)];  % match z at p2?

Almost

Let P1 = (x1,y1), P2 = (x2,y2) are two points on the line. At P1, the difference of z between plane left and plane right is:

0 = deltaz@P1 = (a1*x1 + b1*y1 + c1) -  (a2*x1 + b2*y1 + c2) 
0   = [x1, y1, 1, -x1, -y1, -1]*[a1 b1 c1 a2 b2 c2].'  

The same calculation at P2 leads to
0   = [x2, y2, 1, -x2, -y2, -1]*[a1 b1 c1 a2 b2 c2].'  

So if Aeq is defined as
[x1, y1, 1, -x1, -y1, -1;
 x2, y2, 1, -x2, -y2, -1]

The constraints is
Aeq*sol = 0

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <huret0$iph$1@fred.mathworks.com>...
> Let P1 = (x1,y1), P2 = (x2,y2) are two points on the line. At P1, the difference of z between plane left and plane right is:
> 
> 0 = deltaz@P1 = (a1*x1 + b1*y1 + c1) -  (a2*x1 + b2*y1 + c2) 
> 0   = [x1, y1, 1, -x1, -y1, -1]*[a1 b1 c1 a2 b2 c2].'  
> 
> The same calculation at P2 leads to
> 0   = [x2, y2, 1, -x2, -y2, -1]*[a1 b1 c1 a2 b2 c2].'  
> 
> So if Aeq is defined as
> [x1, y1, 1, -x1, -y1, -1;
>  x2, y2, 1, -x2, -y2, -1]
> 
> The constraints is
> Aeq*sol = 0
> 
> Bruno

Thanks, I think I understand more how to program equations sytem in Matlab.
Tell me if I'm wrong, but I just understood that the part of the code with Aeq is just there to set a common separation line between the two planes. When I first read your code, I tought the intersection was set to the line defined by the user, but I now think the line traced by the user and the separation line just share X an Y coordinate but have not the Z. Right? 
If what I just said is right, the code is even better than I tought for field nivelling...

I have another question, if don't bothers you:

A future fonction of my nivelling field GUI is that the user would be able to set a cut/fill ratio of the field. This cut/fill ratio means that the planes must be 'lower' than the actual regression make them because in a field, when you cut some soil from a higher part than the plane, it takes actually a bigger amount of soil to fill a hole of the same volume than the higher part. The standard cut/fill ratio is 1.25. 
Do you have an idea of how can I incorporate this cut/fill ratio into my regression?  Tell me if you want more details.