function demo()
%This function 
  p=0.6; n=20;
  choice=0;
  [pts,xMin,xMax,yMin,yMax]=random2pts(p,n);
  im=printInit; im=printPts(im,pts);
  figure(1); imshow(im);
  while(choice~=10)
    choice=menu(['n=',num2str(n),' p=',num2str(p),' To resume simulation, press a button'],'Generating new spatial points','Show true bounding box','MAX-algorithm','BTC-algorithm','LR-algorithm','DMX-algorithm','BHM-algorithm','Modify the ratio of outliers','Modify number of spatial points','End simulation');
switch(choice)
      case 1
        [pts,xMin,xMax,yMin,yMax]=random2pts(p,n);
        im=printInit; im=printPts(im,pts);
        figure(1); imshow(im);
      case 2
        im=printBB(im,xMin,xMax,yMin,yMax,0,0,0);
        figure(1); imshow(im);
      case 3
        [xmMax,xMMax,ymMax,yMMax]=MAX2BB(pts);
        im=printBB(im,xmMax,xMMax,ymMax,yMMax,1,1,0);
        figure(1); imshow(im);
      case 4
        [xmMax,xMMax,ymMax,yMMax]=BTC2BB(pts);
        im=printBB(im,xmMax,xMMax,ymMax,yMMax,1,0,1);
        figure(1); imshow(im);
      case 5
        [xmMax,xMMax,ymMax,yMMax]=LR2BB(pts);
        im=printBB(im,xmMax,xMMax,ymMax,yMMax,0,1,1);
        figure(1); imshow(im);
      case 6
        [xmMax,xMMax,ymMax,yMMax]=DMX2BB(pts);
        im=printBB(im,xmMax,xMMax,ymMax,yMMax,1,0,0);
        figure(1); imshow(im);
      case 7
        [xmMax,xMMax,ymMax,yMMax]=BHM2BB(pts);
        im=printBB(im,xmMax,xMMax,ymMax,yMMax,0,1,0);
        figure(1); imshow(im);
      case 8
        p=input(['old ratio p=',num2str(p),' and new ratio p=']);
      case 9
        n=input(['old number of spatial points n=',num2str(n),' and new number n=']);

    end;  
  end;

function [pts,xm,xM,ym,yM]=random2pts(p,n)
%This function select randomly n spatial points, among these a proportion 
%of p are inliers and a proportion of 1-p are outliers.
%pts is a matrix and the first row lists the x-coordinates of the spatial 
%points and the second row lists the y-coordinates of the spatial points.
  tailleIm=[120;160];
  xm=5+round((tailleIm(2)-15)/2*rand(1));
  xM=tailleIm(2)-5-round((tailleIm(2)-15)/2*rand(1));
  ym=5+round((tailleIm(1)-15)/2*rand(1));
  yM=tailleIm(1)-5-round((tailleIm(1)-15)/2*rand(1));
  tirU=(rand(1,n)>p);
  ptsU=[rand(1,n)*(tailleIm(2)-1)+1; rand(1,n)*(tailleIm(1)-1)+1];
  ptsF=[rand(1,n)*(xM-xm)+xm; rand(1,n)*(yM-ym)+ym ];
  pts=ptsU.*[tirU;tirU]+ptsF.*[1-tirU;1-tirU];


function im=printInit()
%This function creates a grey image
  tailleIm=[120;160];
  im=0.5*ones(tailleIm(1),tailleIm(2),3);

function im=printPts(im,pts)
%This function prints the spatial points with white crosses
  satureInt=@(A,N)max(1,min(round(A),N));
  if (max(max(max(im)))<=2) val=1; else val=255; end;
  [M,N]=size(im);
  pts=reshape(pts,2,prod(size(pts))/2);
  for k=1:size(pts,2)
    im(satureInt(pts(2,k)-2,M):satureInt(pts(2,k)+2,M),satureInt(pts(1,k),N),:)=val;
    im(satureInt(pts(2,k),M),satureInt(pts(1,k)-2,N):satureInt(pts(1,k)+2,N),:)=val;
  end;

function im=printBB(im,xm,xM,ym,yM,cR,cG,cB)
%This function shows a bounding box whose extreme points are of
%coordinates (xm,ym) and (xM,yM) of colour cR,cG,cB
  tailleIm=[120;160];
  xm=max(1,round(xm));
  xM=min(tailleIm(2),round(xM));
  ym=max(1,round(ym));
  yM=min(tailleIm(1),round(yM));
  im([ym yM],xm:xM,1)=cR;
  im(ym:yM,[xm xM],1)=cR;
  im([ym yM],xm:xM,2)=cG;
  im(ym:yM,[xm xM],2)=cG;
  im([ym yM],xm:xM,3)=cB;
  im(ym:yM,[xm xM],3)=cB;


function [Xmin,Xmax,Ymin,Ymax]=MAX2BB(pts)
%MAX-algorithm
  tailleIm=[120;160];
  Xmin=min(pts(1,:));
  Xmax=max(pts(1,:));
  Ymin=min(pts(2,:));
  Ymax=max(pts(2,:));

function [Xmin,Xmax,Ymin,Ymax]=BTC2BB(pts)
%BTC-algorithm
  Xm=mean(pts(1,:));
  Xs=std(pts(1,:));
  n=size(pts,2);
  nX=sum(pts(1,:)>=Xm);
  Ym=mean(pts(2,:));
  Ys=std(pts(2,:));
  nY=sum(pts(1,:)>=Ym);
  Xmin=Xm-sqrt(nX/(n-nX))*Xs;
  Xmax=Xm+sqrt((n-nX)/nX)*Xs;
  Ymin=Ym-sqrt(nY/(n-nY))*Ys;
  Ymax=Ym+sqrt(nY/(n-nY))*Ys;

function [xMin,xMax,yMin,yMax]=LR2BB(pts)
%LR-algorithm, LR for learning. 
  tailleIm=[120;160];
  n=size(pts,2);
  sature=@(x)min(1,max(0,x));
  fun=@(x1,x2,x3,x4,x5)[1 x1 x2 x3 x4 x5 x1.^2 x2.^2 x3.^2 x4.^2 x5.^2 x1.*x2 x1.*x3 x1.*x4 x1.*x5 x2.*x3 x2.*x4 x2.*x5 x3.*x4 x3.*x5 x4.*x5 x1.^3 x2.^3 x3.^3 x4.^3 x5.^3 x1.^2.*x2 x1.*x2.^2 x1.^2.*x3 x1.*x3.^2];
  [x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6]=attributs(pts,tailleIm);
  expliXm1=fun(x1,x2,x3,x4,x5);
  expliXM1=fun(x1,x2,x3,x4,x6);
  expliYm1=fun(y1,y2,y3,y4,y5);
  expliYM1=fun(y1,y2,y3,y4,y6);
  if (n<15)
    [xMin,xMax,yMin,yMax]=MAX2BB(pts);
  else 
    nEq=min(35,n);
    [PXm,PXM,PYm,PYM]=donnees(nEq);
    xMin=tailleIm(2)*sature(expliXm1*PXm);
    xMax=tailleIm(2)*sature(expliXM1*PXM);
    yMin=tailleIm(1)*sature(expliYm1*PYm);
    yMax=tailleIm(1)*sature(expliYM1*PYM);
  end;

function [x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6]=attributs(pts,tailleIm)
  x1=median(pts(1,:))/tailleIm(2);
  y1=median(pts(2,:))/tailleIm(1);
  x2=mean((pts(1,:)-x1).^2)/tailleIm(2)^2;
  y2=mean((pts(2,:)-y1).^2)/tailleIm(1)^2;
  x3=mean((pts(1,:)-x1).^3)/tailleIm(2)^3;
  y3=mean((pts(2,:)-y1).^3)/tailleIm(1)^3; 
  x4=mean(abs(pts(1,:)-x1))/tailleIm(2);
  y4=mean(abs(pts(2,:)-y1))/tailleIm(2);
  data=sort(pts(1,:)); x5=data(5)/tailleIm(2); x6=data(end-5)/tailleIm(2);
  data=sort(pts(2,:)); y5=data(5)/tailleIm(1); y6=data(end-5)/tailleIm(1);


function [PXm,PXM,PYm,PYM]=donnees(n)
%calculé avec cptMax=50000;
  switch(n)
    case 15,
      PXm=transpose([0.3061101573    3.3165543713   21.0545444128   -9.7488725728  -11.4424159993    2.3591327500    6.9800063044   98.3377357635  -25.4778005964   68.3671382904   16.4828513028    2.2410080947   14.0410840360  -26.2440978355   -3.9019982851  -43.9965246474 -211.2328498128   59.2178755856  100.1383386201  -38.4561137583  -35.9281120481  -12.8632296295    1.3766313775  -10.0344686324   23.3005125647   -9.3706676935   98.7313045721  -71.3125918400  -80.6674394450   76.0400516513]);
      PXM=transpose([-0.2228974261    0.8175066750  -12.1637014865    6.2688997037    8.6500528320   -0.8441697309  -21.6324132390 -242.7012621239  -61.5035997653  -82.3461367695   -1.9255792451    3.7455041469   -6.0294311476   27.4503578517    3.5712432770  273.9604152380  266.1669338504  -42.3796116850 -159.1970199041   24.0478443392   24.2222854149   21.3966875951  -12.3947796404   49.2207372075   -6.0445891338   -2.0806012618 -113.2400139132   93.2521095196   96.7004642995 -131.3918973119]);
      PYm=transpose([3.6300430404     3.8614408256   108.2994089924   -60.2228792826   -80.0849850464     4.5211275069    -2.2093522980   985.4562653924   310.5521027402   468.4258560112     6.0337129560    14.9910133847     4.6337621736   -23.7816781933    -1.1293118354 -1118.7377516806 -1348.3224598809    89.0332961894   758.3698345337   -49.9365648300   -61.9963449939    -3.8695315801    -5.1622649778   -18.8509371663    38.9005569801    -4.7165628932    73.5623017992   -76.2098946621   -71.2644378933    96.2826684570]);
      PYM=transpose([1.0932232754  -0.1487187227   3.8017502166  -1.9456433957  -2.0258742594  -0.4700399247   0.5581014621 -10.1472539050 -29.6930073152   5.5426209767   0.6068515625   3.2727539744  -2.8985172052  -2.1270184815   0.3494986777  45.0327448762 -13.1928504472  -5.1723769832  -0.9068442945   3.9464870093   2.2924145122  -0.5556705618  -3.6677508592   2.1077414495   5.0628019845  -0.4901830146   3.3474077189  -4.3311828169  -3.4211013000   6.8461359354]);
    case 16
      PXm=transpose([0.3104935162    3.0655658636   20.7419156106   -9.4682825966  -11.3870512530    2.6176433552    8.1262130438   87.4420819525  -42.2076625740   68.7290865530   17.0352481911    7.2966149225    5.9135001521  -26.4210439396   -4.8483278676  -15.1959297797 -211.9508442616   65.7486270226  100.3390943621  -42.4234404546  -37.9318329075  -12.0642317987   -8.9020372384    2.3020490412   22.8439394847   -9.8295516021   82.5591110820  -54.3117583707  -66.1214197698   60.8864458942]);
      PXM=transpose([-0.1613179833    0.7581736552   -7.7519316083    3.5280114209    7.2218256128   -0.8942654441  -25.1726180253 -164.7030267391  -25.5248311056  -69.0780866909    0.5713043462   -0.5174203929   -1.8197160886   32.4289580396    5.0066627055  171.4632150444  194.2830722588  -26.0867549001 -115.1802856973   14.3077951269   14.9351516923   20.4667130923    0.0158341281   35.3553027153   -0.7664184714   -3.6384766729  -95.4066852752   58.5702346059   87.2637495519 -106.6932087424]);
      PYm=transpose([3.9787760503     4.3372973860   120.5776151470   -67.7712164139   -87.9620241908     4.4667729580    -0.8786734654  1093.5240883168   339.5171767004   508.5695739487     6.9473934646    -8.8574116092    28.1799224999   -23.4433412300    -1.0655456143 -1238.0068548214 -1463.8451795525    90.2328448158   818.7641692971   -49.2040275951   -63.6890205017    -9.8820610832    10.8160950259   -37.9448896280    38.5674447460    -5.8863564807   126.5143069499  -124.8628948144  -113.2823683302   141.9183290158]);
      PYM=transpose([0.8096467414   0.0668721869  -4.2160249478   3.4950068091   2.8945428081  -0.5455328207  -0.1170522077 -65.2114869144 -62.2868788785 -10.4513454886   0.7675978181  16.8238146959 -10.9477163233  -6.5069358274   0.1134556480 126.0487916733  36.9282814555  -4.5122720992 -36.4803684155   3.9215494052   2.0541589842  -0.1870124615   2.4755401687  -6.4735252627   9.4474548303  -0.5989751627   6.6524193166 -18.7322140438  -7.6475295991  23.5384195887]);
    case 17
        PXm=transpose([0.2789511051    2.7926492324   19.2574452181   -9.2803976385  -10.2904050201    2.5145833748    6.8058935977   65.1098265152  -38.7269341687   60.6159124915   16.8176668618    1.2851008156   10.2706132439  -22.4060300967   -4.0507953937   -2.2217762056 -182.9237715107   62.0700207794   86.7074345452  -39.1872041754  -36.9373213516  -12.3341708038   -1.2255971829   -2.0272854612   20.9430993383   -9.9417844750   87.4010401477  -59.8584727877  -67.8609631592   58.1451370802]);
      PXM=transpose([-0.1863732489    0.9225630672  -12.6606668676    7.3690437318    8.6874398185   -1.5980822748  -26.7319547440 -212.8675286755  -33.6580399523  -76.0238994177    2.4789013332   -4.5930736977    3.3905787661   35.0481638661    4.4722360053  224.4066884526  237.8347620514  -27.5875136116 -147.9271382234   15.4974038315   14.5987950054   27.3581376905  -14.1350701360   53.5195859072   -2.7290568318   -4.3710805092 -128.6624345959  107.5730806453  100.3212695766 -145.2883415945]);
      PYm=transpose([4.6429833012     3.0218737640   145.4843128258   -78.9726428933  -102.5219883830     5.8182430132    -1.0213302600  1288.4798157838   407.3109440257   591.8845836176     8.7516133059    12.2445559652     4.8422530850   -20.4871001925    -1.1678341792 -1452.8006166655 -1754.7299741055   109.2776540652   973.8328328250   -59.1128751798   -79.0652136829    -4.4548484044   -19.2574807556     1.9799976079    67.5918343190    -7.0902845405    62.5501319516   -70.0272029715   -54.6430466913    74.2205434318]);
      PYM=transpose([0.6804464768    0.1676708122   -5.0968667869    4.0863249801    4.4143289248   -0.4943141647   -0.0923315453 -123.3565289012  -88.7620248110  -28.2420298525    0.5519200037   17.1924996374  -11.0988839111   -7.2105474467    0.0761468813  208.3912976921   95.4796997514   -8.3638130802  -79.6978438798    6.7024032128    4.0731209164    0.3443710916   -2.7730669006   -3.7993601911   13.4767102061   -0.5333872751    2.8505611741  -10.7467334136   -6.0458924128   17.7946611098]);
    case 18
      PXm=transpose([0.2915043796    2.1230544934   16.9117876681   -7.5332696289   -9.8345127463    2.6742517606    5.1016327149   29.8449341498  -72.6093663140   56.0940920044   20.0458837514  -28.3668757669   30.1616808964   -9.7712853351   -6.0262125098   67.7971783956 -150.6881077750   79.1770881890   59.3668865704  -48.5061600177  -43.6124814421   -8.2630152363  -17.6896296664   22.3731455578   14.1300238908  -13.0090962900   77.4517352528  -22.8840261090  -71.8965808105   31.7589702306]);
      PXM=transpose([-0.0654083395   0.6078474260  -5.5742593935   2.2313758060   6.3454551878  -1.7616315320 -28.9036854813 -73.0313029656  22.8355308606 -56.7551930708   4.2203290403  -3.4905601274  -2.8287085965  39.1285440306   6.5091634670  40.5034897254 126.7382523245 -21.9423639426 -69.7100166417  14.9164753009   8.4490299683  16.4719207726  38.8784346669  12.1236369058   0.7599999902  -5.5227725202 -58.8157577265  -9.5574229057  76.3337070214 -77.7075634757]);
      PYm=transpose([4.6546201456     4.4423438051   150.1096264151   -82.1565495320  -105.7449929806     6.0887549252     1.5878019228  1286.0926055396   372.5644877889   610.6518969898     8.7791793597     3.8568029690    17.3745066150   -31.7330092672    -1.0628421167 -1408.4142773396 -1780.7617720224   115.4488443152   975.6282866962   -63.5107738827   -82.2434458898   -10.6444936340   -14.5689226516    -4.9161691751    72.9985622208    -7.2970242860   108.4821429247  -108.3788981236   -89.5968572625   109.7243251872]);
      PYM=transpose([0.9987142567   0.1450423685   0.8711407945   0.9644327229  -0.1107976036  -0.6869709210  -0.8325775193 -77.4005048824 -82.7340178643  -7.5714192108   0.7721719677  21.4891758768 -14.3339564188  -6.7994276613   0.0814587696 164.7927458750  29.9786824901  -9.9311994213 -43.5232267481   7.9556014304   4.8178112276  -0.3086432567   4.0938759496  -9.8733045796  15.5107145990  -0.6399606656  10.0259967121 -28.2791379440  -9.3681017859  31.5127576770]);
    case 19
      PXm=transpose([0.2611682701    2.4149126555   16.0373090383   -7.2532140305   -9.2264771426    2.4597689226    5.7520153004   13.7736080519  -73.6074409589   53.0080689853   17.0966657392   -4.5194731261   11.3544272104  -17.2834421602   -4.7813116312   80.6454815850 -143.4915946050   65.1748081502   60.1784128177  -40.7413522473  -36.6535584551   -7.8458163049  -18.6717728683   19.4936133582   15.8755762441  -10.6987075511   62.2458267400  -18.7979355333  -55.9327779701   27.1094731780]);
      PXM=transpose([-0.0832966329    1.2162276328   -5.3640534263    2.9683152251    6.0467292151   -1.6334069144  -29.4106991305 -110.4404412373  -11.5165328772  -59.9277472682    3.7956116765  -19.0068418087   19.8217370088   39.9307794663    4.9071521670  116.2564479142  150.1877528503  -17.6673769865  -98.5784796809    8.5921065255    8.7878634566   21.2219787993   25.7972339714   13.2655348033    2.6486809131   -4.5959844780  -59.9975119079    6.3034902764   53.6460245013  -65.1323893709]);
      PYm=transpose([4.4478541322     5.0446608302   143.3425497410   -76.8396333975  -101.5733878576     5.6710854009    -0.3661339412  1230.5345458031   351.7462342244   589.6476151600     8.4128600979    19.6529686752     9.9175607151   -37.2098757706    -0.3139138328 -1336.0160426785 -1719.8790715646   110.1995438039   928.6448747870   -60.7496291887   -78.5507904533    -8.6149251816    -6.4315687084   -18.1285399690    80.9865245106    -7.3215807683   104.3633681777  -120.7395136345   -88.3575736602   123.0693554183]);
      PYM=transpose([0.6513382316    0.0473758965   -6.8546254860    6.0310822327    5.5306256086   -0.4609629873   -0.5891202663 -197.5737392584 -166.6247517606  -45.8643037485   -0.1245347365   28.4279994619  -20.8786034243   -9.8527753616    0.6118414822  362.5038864095  154.8553105746  -21.7171097143 -136.5411502090   17.1484561279    9.9634483896   -0.7684515288    0.2603891719  -10.2813570238   21.3973039432   -0.2990307831   13.2979891260  -36.2266910575  -12.0298749459   43.3829343776]);
    case 20
      PXm=transpose([0.2361013269    1.7007433089   13.1318176453   -5.2852360701   -8.1194451370    2.5793343392    6.0474653222    0.8346914814  -78.4772251789   48.7848206041   18.1637866097  -18.6817494262   20.5616141851  -11.4413193769   -5.3655386097   98.5259594202 -121.7752610757   75.2832715238   43.9093316033  -46.4955516256  -40.1629876311   -6.0787229900  -27.3579300893   29.5475375498   12.3505716130  -12.2030335903   50.2161818400    3.2897740608  -51.1000811019   11.6747847115]);
      PXM=transpose([-0.0371166564   0.9473126407  -1.0920611385  -0.1079470019   4.7821454184  -1.4061945608 -36.6175041838  -9.3695297592  32.7560396253 -49.2149605640   4.7998063274  -5.5889172046   1.1123768002  49.9200376626   5.9790980265 -16.0053589344  64.8841003276  -7.7078981320 -37.0533011660   3.8741476906   2.0849608395  21.3185708865  56.9473599267   0.1991432302   6.4105031615  -5.1164112585 -50.9940810141 -55.8503503321  72.0875255489 -50.4104818399]);
      PYm=transpose([4.0727190755     5.0402461847   145.6009907763   -77.3287043463   -99.1950941152     6.9878086820    -2.2688635693  1239.1265091880   318.2458426059   587.2553320839     9.5370248380    17.3342749146    14.4125387606   -33.3627067926    -1.2781404104 -1291.5146885904 -1754.5893859319   126.1540303721   924.6609065515   -68.8279725364   -90.3010066832    -6.2426870481   -27.6024903801     0.9502249056   107.7290225235    -8.1855077845   100.6485803332  -116.7679920038   -90.6526343729   121.2675758002]);
      PYM=transpose([0.4885602032    0.1295107813  -13.1724158851    9.5862216430    8.9701394879   -0.5963258090    0.8983193769 -231.1794848854 -163.5935379861  -58.2895957879    0.4021261149   29.2121913646  -22.7229389967  -11.5357541590    0.0487403438  381.8803791083  202.4333465605  -17.4797963830 -157.4068764585   14.7627794209    8.0176247724   -2.6167809020    7.0808045980  -13.2435246075   16.2127126408   -0.5591236559   15.1845145252  -38.1608420014   -8.1958495638   38.5601404410]);
    case 21
      PXm=transpose([0.2641847015    1.6285747821   15.4979798827   -7.1875397550   -8.8066817588    2.8237172558    3.7615018139    8.1686625474  -92.0544505393   48.6913246813   18.1131278764  -22.7747218365   19.4641729685   -5.7964151018   -6.3472429249  109.4341065784 -133.5837654063   77.9332481481   52.7135141592  -47.5778083102  -40.5183972574   -6.8389088495  -27.0581957995   24.5388683889   11.5799815010  -12.3177141584   67.9454120857  -15.4748563956  -64.0746688145   31.8188658427]);
      PXM=transpose([-0.0745090354   -0.3299712857   -7.1223613118    2.8531431311    7.1321604148   -1.4813042721  -35.0156990794  -44.1679182606   65.2230087474  -61.3047653532    5.2608281761  -48.0545172105   35.5864551572   61.2391844202    5.5880317524  -26.2517398550  133.9058750583    2.2301690904  -74.0692351736   -4.1100099370   -1.0784978943   24.8196598172   47.3376980988   22.8955878008   -2.7837900326   -5.1603725904  -76.6661556566    3.5929365034   83.2747857750 -109.4467604396]);
      PYm=transpose([5.0575117299     4.4431802104   170.5786454127   -92.2558567963  -116.8417499649     7.2868669189    -0.6790713561  1480.5014648954   435.3557177937   667.1230735891    10.9073940370   -13.7183569555    37.0521252898   -23.7892772859    -1.3817911806 -1630.1307683188 -1997.9185519439   122.9244495456  1085.3487452931   -64.6799208142   -91.8098525251   -11.5380557428     3.6260461944   -24.6673140968    94.0802333743    -9.2001307301   137.6884798537  -158.6444033582  -111.4892368422   149.1652644052]);
      PYM=transpose([0.1267644597    0.0447607225   -7.5548903784    5.1944278009   10.7033756999   -0.4398463702   -0.2639398471 -168.2174090111 -143.8922660823  -55.2379846504   -0.0983572969   22.1993018643  -12.3420688354   -8.9466209640    0.0125507886  312.0348623951  136.4641705304  -15.8530690676 -118.9933701653   12.0627087762    8.1434928345   -1.6991434277    1.2828052250  -11.7177656540   31.4045289245   -0.1532273068   21.8432584145  -44.0421524829  -19.1389167599   47.2428111777]);
    case 22
      PXm=transpose([0.2470224462    1.3900153590   13.4211348663   -5.5907036198   -8.1945055545    2.9339482885    5.6497734415   -1.8529508326  -93.6966403043   49.1855981743   17.0843125715  -14.5545635920   15.4451235477  -10.4383049968   -5.5314054262  121.2718630035 -126.7594509967   76.0092521334   44.7376247379  -45.4298381194  -40.0713099192   -7.2965103550  -26.3624559226   21.8568053760   11.8646239806  -11.8389850976   62.5225338234  -13.5609614011  -61.3095373983   32.6043443323]);
      PXM=transpose([-0.0248267243    0.3329912976   -0.2524272773   -0.9936404066    4.5101922633   -1.0472591929  -38.1615303154   43.8881399714   79.5037479809  -46.8474852008    4.7403485474  -26.1069549743   15.1227503818   60.2261552121    5.8877806985 -112.8409999143   43.0631734849    4.5150771836  -16.1518221490   -6.0335734518   -3.2507119742   19.7654793352   78.2101335588   -9.1717357081    4.0253418706   -4.3948790136  -38.5831654196  -82.0428810677   73.1989502165  -50.2228739785]);
      PYm=transpose([4.2209588426     5.1091142700   152.8437270287   -83.1350393388  -102.6643709766     6.7673807818    -2.6085013802  1340.5072761759   372.8658557378   603.7477774583    10.3003583287    -5.6844631951    39.2236257991   -28.8222835195    -0.9535700693 -1445.0052868278 -1822.6867356988   118.8896537408   973.8419449590   -63.1332736014   -87.6923496096    -9.7316166815     2.3456039219   -20.4131736705   101.6452494922    -8.9711723253   139.3545929612  -165.7670513980  -117.6200464001   152.5639061874]);
      PYM=transpose([0.4535733260    0.0157058130   -6.7398903426    5.4880071229    8.7255446756   -1.0733272598   -0.8204739160 -217.8195395003 -194.4091409611  -63.8102459953    0.3464696770   20.4581880820  -11.9960323550   -7.0911788006    0.1764300755  410.9318234169  181.9753664961  -25.2508892604 -160.5614326570   19.5473210402   12.6915219860   -2.2744822465    6.9801708169  -18.2544018998   29.7052210419   -0.4715962407   31.7843295654  -61.2872130464  -26.7398072713   65.8834209266]);
    case 23
      PXm=transpose([0.2182621726    0.7309031431   10.6641334521   -3.9781897424   -6.9463661455    2.8780696804    2.0681316077  -10.9439176708 -111.1970943457   40.8011489219   17.0255448888  -25.8271492855   17.3987318194    1.6842620136   -6.3977464572  144.6291328825 -107.8944977513   81.7327991583   37.8414933167  -50.2686674631  -40.1666735806   -5.9012897395  -13.4996870309   12.7566175068    6.9730877865  -12.1838544448   66.9184131514  -32.0396039704  -60.6750208018   45.8402262995]);
      PXM=transpose([-0.0024159925   -0.4539494827   -0.4733606092   -1.8647217111    4.7609209124   -0.9301176774  -38.6923414162   70.8039551222  105.2779138475  -48.3381776004    4.9874086653  -55.3319223347   37.6380458544   68.9794752321    6.5963455879 -168.7053529981   49.8440085237    9.8577891491  -14.8989601216   -9.3342974734   -6.4313908399   19.6385711747   93.8069035949  -16.4605803465   -0.1895496985   -4.4754241606  -25.5423481128  -98.3962190903   64.3570246398  -45.5021567710]);
      PYm=transpose([4.7163993297     4.8367199083   163.0650626808   -90.0633916415  -110.1109204571     6.1463526384    -3.8057458684  1482.3003575677   415.5138941374   650.1523960028     9.7600540442     7.1939336548    26.1362865296   -27.7678264823    -1.5846525982 -1615.1888448380 -1990.8841356935   116.3285793845  1091.4881080996   -63.2673486673   -82.0827643603   -12.8678981256    22.9114419561   -33.1355927750    93.7205880087    -8.8620165176   160.9344181243  -212.6730128273  -122.6690991864   181.8610164771]);
      PYM=transpose([-0.8385037600   -2.0348751229   -5.7093958761    7.7608174910   14.5424804634    2.1813625044    0.7333404476 -206.7412124096 -261.6333390024  -76.1232339176   -4.8313237554   18.1765859578  -23.5447334424    2.6749325038    0.8060471716  499.4237688594  147.1952206252  -31.9576229586 -170.0201796906   23.4836999632   15.5892756232    8.3775784134  -77.3950118323   59.2402382036   51.4610600881    2.1178311461  -63.2177572806   69.5119345980   32.9483593673  -28.9175815328]);
    case 24
      PXm=transpose([0.2111812853    0.8378079931    9.8610708027   -3.2419194611   -6.7624590275    2.8037667248    2.6621372981  -30.6972439143 -113.2530634991   37.8879187369   17.2017059262  -36.3619380633   29.9552755394    1.3348308170   -6.1068515724  165.6785076689  -83.0639562035   79.2155779140   16.2560994308  -47.7997231300  -39.5096579308   -2.9859895111  -31.0170618352   30.7210362932    5.5947247523  -12.7255835273   52.0605348108    6.3723636695  -61.4018703576   20.8396750407]);
      PXM=transpose([0.0405862406   -1.0038986630    3.5404569687   -4.8946740448    3.7781464616   -0.6881111540  -43.9019461984  165.4041312760  149.8125823564  -44.1460531428    4.5553036999  -49.5331281779   23.9668879295   78.8997780474    7.6019333418 -300.9417650304   -9.6655951567    7.8185785280   35.1271869556   -7.5141947628   -7.4004031066   17.8545218904  131.5027109456  -41.1076561857    0.1994647386   -3.9875650616   -7.1149193527 -166.9402684737   73.8220391777  -16.7999461580]);
      PYm=transpose([4.7234970679     4.4037522065   176.7770033556   -98.3027227333  -114.7812755960     5.8742383864    -1.4058575408  1642.3674032815   461.7594161875   696.8973786852     9.9163004363     5.8644616122    27.4453553383   -30.0480708504    -1.4803669548 -1793.5067360322 -2183.4506502557   121.3915641822  1204.7865512199   -67.6661096137   -82.8296149555   -13.2304238606    11.0263077044   -30.1138438153   107.2541653701    -9.5349038293   156.0254374276  -206.1682435742  -123.0494351270   187.5973837938]);
      PYM=transpose([0.1730220223    0.3220559856  -14.2649794339   11.2211438439   11.0782296347    0.1124550601   -1.4170355664 -281.9790490483 -211.8588205751  -78.8031637731   -1.3784511201   34.8105916544  -23.9714853166  -11.1159673094    0.2069712559  484.9863821239  250.0714882027  -26.1887764193 -204.1322670836   19.6033923442   12.9850334257    1.2855399625   -4.8517263970   -8.1438075434   27.0114390129    0.4748045815    4.6806529197  -29.9441676850   -9.7286223491   43.6466378013]);
    case 25
      PXm=transpose([0.2323071735    0.6977700777   11.6103844630   -5.3170451336   -6.8921145773    2.4715229226    3.9937553929    6.2759230741  -84.4920557102   40.3141030804   16.4566138495  -27.6158880304   20.4985597715   -1.1112477455   -5.6012228803  100.4852271903 -107.1405039781   73.4830554649   40.7719275641  -44.0429868102  -36.4925001479   -7.3533893408  -23.7264534279   19.6938662123    5.4029053332  -12.5096186425   64.1871076313  -12.5097608543  -61.2408950681   32.3561611458]);
      PXM=transpose([0.0323625806   -0.6506887745    3.3862418437   -4.2102239529    3.2874999528   -0.4390565586  -42.0498486567  144.2710914764  121.3091293754  -44.4589057564    3.7338967253  -71.4840316524   45.0402274684   78.6504020959    7.4623854805 -244.8992663813   10.4551880901    9.3117543923   11.1682838654  -10.2986335484   -7.2193338134   10.4133709935  154.6751531671  -51.4409508395   -0.9783916495   -3.2812015518   47.3285597738 -224.6495576892   42.2759985277   11.6799408571]);
      PYm=transpose([4.4126527774     4.2587381314   171.8613700131   -97.5595460857  -109.1475989493     5.0651580727    -6.1499093719  1740.8843242843   537.0457942382   677.4906302906    10.2187497017     8.8571877345    24.7554451305   -20.9639594668    -1.0544654795 -1974.1963396291 -2205.8802247775    94.6304904222  1249.7095294008   -49.3717176559   -70.0275710473    -6.8998848767    11.7642142618   -38.4075235899   101.3490544879    -9.3655565862   134.5162425506  -195.5771284164  -115.3749047468   191.8339767123]);
      PYM=transpose([0.3086782568   -1.5131991534   -7.3649390559   11.5644182727   16.4445376784   -3.6607503111    1.9003721462 -447.4701465657 -444.7026168471 -131.5019227639    1.8792192205   39.5167156649  -32.9407672250  -10.1720766181    0.3753874730  921.2232713500  352.3156201492  -64.8098018015 -357.1970710744   48.3289141331   33.4369284253    0.6051690608  -53.1052102907   37.0371785717   76.0813053667   -1.3743438815  -13.0487769269  -10.9787616705    6.1319640914   29.2147636463]);
    case 26
      PXm=transpose([0.2079618324   0.4034370219   9.1924462837  -3.9689950643  -5.8532004507   2.1781284615   1.8523159586  11.0715616614 -81.6282193532  34.3414933724  14.0626736183 -20.2031552838  13.1054572558   1.5539727042  -5.4711104744  88.4954518183 -95.0082509553  60.9088641955  38.6869194581 -36.9438880097 -29.2934200732  -4.2911298746 -16.7502706876  20.6715696530   1.2531359699 -10.7896343476  49.9801285838 -10.6357670658 -49.9576441045  26.9101255253]);
      PXM=transpose([0.0117350503    0.0692117083    4.8632774762   -4.7406616877    2.5851852855   -0.2253965132  -45.0824664557  150.7334364438  126.7071472656  -38.9542082465    3.9085513105  -48.0588288188   26.7742642615   76.2126006850    6.9737930135 -248.9360027295  -20.8979091790   15.8363960776   29.3015707500  -15.1982532689  -10.0671456284   15.6383078090  141.9581869107  -41.9899851424    4.6525597621   -3.0629872283   18.2989863262 -205.6960364011   63.7743804885   -6.5462762707]);
      PYm=transpose([4.5355354630     4.4091915510   180.9353749576  -101.7891122354  -114.0948787939     6.4378199850    -4.6968754323  1744.8119223504   514.7877512869   705.0299262104    10.4156182120    20.4404461063    18.9303886452   -29.2757206109    -1.0885881515 -1938.2380346889 -2275.3311525684   111.4575349906  1268.8195696730   -59.4345352958   -81.8441772480    -5.6743896440    -8.6340462713   -12.6436923384   121.7065843401    -9.5083082412   116.3002977295  -172.0808826165  -101.2170437474   163.2131334139]);
      PYM=transpose([-0.2740637016   -1.5792566557   -9.3006257536   11.6817314465   14.4583012373    0.2870401354    2.6064685027 -392.3604480982 -371.8441931827 -121.6251898053   -3.2866635778    3.5721030592  -15.9366290483    3.2363693554    0.8321287460  788.0057126566  346.0350480834  -48.0463709652 -323.3705756762   35.0236514448   24.6172166695   11.9712230885 -116.1529518546  102.1424363541   47.7928212169    1.3207327162  -99.8767629780  139.8155085358   56.1029054648  -87.7747483512]);
    case 27
      PXm=transpose([0.2210517986  -0.1407714003   9.4961347712  -4.6295471063  -5.7425316160   2.2819181228   0.7827892999  23.0771538295 -78.6101173964  33.8018163282  15.7947106298 -35.5313123357  22.2520110897   9.4197303063  -6.5363717003  75.4465009501 -99.1936726180  71.4821049127  43.0316282266 -42.6745068148 -33.4592674489  -5.9945157121  -7.5824144215   6.5762389192  -0.9289549621 -12.5879528704  71.1659164139 -33.1176686106 -65.7339318263  49.2467009758]);
      PXM=transpose([0.0206957174   -0.4022208370    4.9510317660   -5.1566535919    2.9259393695   -0.4126866438  -43.9184211796  169.6551030864  112.5426631388  -43.4959853680    3.6169172478  -79.3433027620   46.6529650095   82.6615251493    7.7595856456 -254.2997615590   -3.4569425771   10.9511918246   18.0929916207  -10.8829739286   -8.7405873649    7.8952262394  173.7163338612  -68.3396879506   -0.0019640995   -2.9699353844   76.3360348670 -277.5849269299   31.5092940553   46.5215432237]);
      PYm=transpose([4.7636784848     3.9047855182   174.4387685407   -98.5122898227  -112.9498516595     5.3851323958    -4.0463181700  1759.3064322675   549.8479341090   694.5718543510     9.6986744716    25.0150618306     7.6988502728   -26.6980229871    -0.2688269581 -2000.4993134896 -2241.3447225404    95.3488775855  1273.4901792588   -50.3568430333   -71.2521079558    -8.5202784401    11.1463031480   -26.6722727994    99.4603478872    -8.9337069537   117.3939125524  -188.3892721469   -87.6192369573   165.3483155494]);
      PYM=transpose([0.4639282637   -0.7229374748   -5.4908131356    7.1942992986    8.8491705883   -0.7746736795    1.6760016795 -423.6429306711 -348.6013956780 -111.4501361467   -1.4931014668   22.9865155128  -17.8153331490   -7.5947821128   -0.0500871067  770.8700477914  363.6080025599  -47.3283206470 -318.6443001891   35.4913754127   24.2557019482   -2.2886567663  -13.0808642579    1.9098815376   40.3384766515    0.5648434904   13.4946784353  -37.8574932981  -10.9076041444   46.6992248441]);
    case 28
      PXm=transpose([0.2253596101   -0.7050613883    9.2003784076   -4.5841083612   -5.5285394967    2.5305881495   -0.2947687720   31.9780209101  -79.0031643050   32.8813276226   15.7465831028  -40.1402890813   22.4343195917   14.9753694302   -7.1677062510   67.5462710220 -103.4992806163   76.8328392359   48.3781425118  -46.3647028410  -34.9478522764   -4.9124870758   -6.2893725285    9.5471484769   -2.8750551672  -12.6993558139   67.7900741062  -31.9842354483  -62.8069202494   45.9524197861]);
      PXM=transpose([0.0370096110   -0.8011916380    6.9498789166   -7.5793719009    2.4007369149   -0.0558466848  -46.9045817918  235.4304845742  182.6735119918  -38.1622940872    3.4685204691  -66.2879484442   32.3713738843   86.9887573114    7.7858054807 -389.6888475182  -50.2650214811   16.7775997993   66.0400544587  -16.0536981453  -11.5408732794   11.7324833657  177.0054040785  -69.0790548312   -0.5319603234   -2.5742158637   48.6862527805 -263.5799342498   58.0306913582   22.7739840186]);
      PYm=transpose([4.5245917024     5.0120554129   165.2445899853   -95.4126929080  -108.4648445781     5.0894943021    -6.5868624122  1771.7912662469   601.9532705357   668.3373444170    10.6177284372     3.0457824167    28.5532343952   -21.1166236067    -0.5841502021 -2088.0199927540 -2165.7288201606    78.2395193830  1256.6170003711   -38.2902855428   -63.5457033852    -9.0033046939    30.9127632052   -29.2337530009    76.3074043509    -9.5522298364   138.8783799441  -199.2616677782  -102.7621139183   158.2890019245]);
      PYM=transpose([0.3791392730   -0.7305112111  -15.4262172515   14.5066945514   14.7538585685   -1.6940891292    1.2734202437 -566.9518777922 -435.0705422411 -157.1604786484   -0.7673918437   28.1642984214  -23.5424194010   -7.6371636498   -0.0847941951  991.7257706009  520.6198397469  -58.6381353238 -436.8458724967   44.7280180642   29.9091130755    2.8475428096  -36.1917999688   23.2022979120   42.7388199563    0.1497633083  -21.3746988952    5.8729051753    8.2367232377   18.7480444468]);
    case 29
      PXm=transpose([0.2218064790  -0.8759342511   8.5741974690  -4.7814040490  -5.0792858327   2.2916958476   0.8077380473  48.3931070511 -48.1790397693  29.6043046985  14.5699525527 -50.0268966331  33.8846631679  14.7949403853  -6.2661089349  16.5051566854 -88.6623356944  65.9142559358  42.0736882346 -38.3462610649 -30.8573356178  -6.6715371133   1.6966810325  -2.2799594432  -5.7547772197 -12.0136889220  78.0046107605 -38.4845495587 -72.3471937613  54.0331854512]);
      PXM=transpose([0.0028666434   -0.7610264542    3.5705568049   -4.7269070906    2.8597192837    0.2343231005  -45.8829407888  192.7019331421  162.9531728256  -41.4498958735    2.8650871599  -73.8667131228   35.3172274624   86.8630764528    8.1721123335 -327.9474095645  -14.0411240334   18.6749284908   40.0646230291  -18.0182392073  -12.5690269666   11.4923265984  167.0274975067  -54.9947927043   -3.6372373514   -2.0711073277   45.4468021482 -245.5119883355   61.1196792562    4.1236326198]);
      PYm=transpose([4.3409010875     2.8062844263   170.6530403168   -97.9231151840  -105.7322894855     4.6878011309    -4.6430948953  1866.0571356507   608.6418344050   665.7837670740     9.8872017254    -5.4082581827    28.3063936045   -10.2246524182    -1.0767446889 -2163.0877035316 -2252.0973753624    68.9867982800  1299.4610048866   -32.1415626811   -56.7971766015    -9.8551112715    11.8858239955   -23.4625582684    96.1913804054    -9.0314830058   134.0302498838  -200.4744190573   -99.6905226371   171.5809766185]);
      PYM=transpose([-0.3405129531   -1.4945096703    1.4850110134    5.7136559591    7.6029159179    1.8846605187    1.4596193460 -391.4396074140 -402.1314669282 -108.4372118274   -5.7426252086   13.2861382628  -23.1849620063    3.6327212436    0.2997313902  825.6222370212  310.7723685222  -57.9538812897 -317.8762308483   42.8295444987   30.0422035118   12.4362332338 -112.4490900837   96.8806890662   52.2724419857    2.3311262412 -101.8037622717  129.5231014925   58.5312488574  -76.2551092315]);
    case 30
      PXm=transpose([0.2055253272  -1.4287075650   6.4487641368  -3.8174370333  -4.1970122617   2.3237323020   0.3026937010  41.4643793053 -42.2306497877  24.6240826422  15.4456233125 -61.6475041697  40.7234352342  20.7811363551  -7.1250590386  12.9980099714 -68.4356454395  69.5694008104  33.4166821042 -40.2168993421 -31.9222683268  -5.3502076411  -0.6220511178  -3.4657697339 -10.1483683328 -13.0346522316  73.9361990627 -24.4461437811 -74.9862902234  50.8526715910]);
      PXM=transpose([0.0231529493   -0.8527778765    7.2905451426   -7.8980374074    1.8343473157    0.5234218132  -46.4687412297  238.0283375808  176.9111887740  -37.7289176230    2.5678341218  -85.6648475073   37.3582411574   90.9328257807    9.1499935486 -382.3372576236  -40.0216610023   21.8550628081   62.4584054516  -19.2395201557  -15.1658706107    4.7653646395  204.3351016970  -76.5138186640   -3.8949914548   -1.9143485957   90.6506517730 -323.0086500871   52.1534727554   36.1063905357]);
      PYm=transpose([4.1065844254     3.6450141632   146.9426918985   -85.1220099171   -96.1575361685     3.6256611614    -8.0056885955  1779.5763859750   636.0364227814   609.6133088372     9.3625031596    -0.2824508703    23.3184026209    -7.7167395795    -0.9906317420 -2151.2239534104 -2057.9019009711    43.7957127354  1223.5965016646   -17.6465214955   -41.0900707406    -5.5517396230    32.9222607128   -29.4112407708    55.8760161717    -8.3974955954   112.1788454692  -181.1083728901   -82.1835734512   147.9821654449]);
      PYM=transpose([-0.4866295129   -1.8867800310   -1.8288690390    3.8663341657   11.8656760146    1.6304492015    2.0427671666 -196.9079670196 -215.0169873001  -77.6487882594   -3.8808240350   14.5255602522  -20.0159375519    1.3447609691    0.4200855570  435.0858935315  158.1911172989  -31.0497890319 -160.4995517674   22.0087448830   15.6746707341    5.2763723852  -67.6056177568   54.7028805163   47.0803944024    1.7365285588  -52.4460513934   60.1425933652   31.9815037844  -33.4820443706]);
    case 31
      PXm=transpose([0.2057541675  -1.5254800354   7.1742378895  -4.2816024017  -4.1295203406   2.1074750516  -2.1338314665  70.5766396954 -39.9745830548  23.2599204156  14.5726425816 -57.4462913896  32.3593934475  25.3505626611  -7.2494294949 -16.3091836955 -86.8503858343  65.7080212833  50.3762115585 -38.7320100555 -28.8297148366  -6.2410647619   6.5493914963  -6.6277855046 -10.3614788741 -12.5934993869  81.5075172624 -46.9528071523 -72.3399854272  63.4780458491]);
      PXM=transpose([-0.0235229519   -1.5193780421    2.8922776143   -5.3884282841    3.7728215819    0.3891310344  -46.0412860586  200.0606117580  159.0973588364  -47.2758938604    2.7290369635 -100.6432138989   46.6761097317   95.1285915802    9.6559592331 -334.9940900784   10.5736787613   18.9110480463   31.3474885351  -16.3565939197  -15.3300410020    6.2436448520  202.0136971601  -79.2954010124   -8.3700100522   -1.8448096695   86.8310303176 -310.1061479283   48.0974304087   37.7038525188]);
      PYm=transpose([3.6864354840     3.9230303694   140.4475634755   -83.9859819646   -89.3996827903     2.8074783196    -7.1516235507  1805.1628945433   632.7406001464   588.8449268400     9.7239094319   -17.5325659432    37.5993553704    -6.2546778524    -0.4114388469 -2175.7947472609 -2034.4214354714    30.2375883070  1221.0866359209    -9.0829852132   -33.7334696722    -5.3828777615    16.5463618615   -26.6688707035    55.4527684844    -8.5293728166   123.0795249896  -178.8560939384  -101.8387227501   168.2684089729]);
      PYM=transpose([-0.8918385652   -2.4888548452  -11.0333922470   12.8821546923   19.8785114486    1.0306097483    3.5334250708 -453.3937528959 -433.9530089373 -150.7742724811   -4.1109080014    1.8409969733  -18.3086175632    9.4308797903   -0.0331925680  913.9291560051  400.9266974742  -50.4600524491 -372.9760934569   37.9989531199   25.5732901103   12.1177847720 -131.9657359347  124.3203416553   60.7650942133    1.7287972983 -117.5483417988  156.4988065219   75.4256502550 -111.0260563701]);
    case 32
      PXm=transpose([0.2092311089  -2.3317280823   6.3430739274  -4.4716543169  -3.5604392279   2.3341962981  -2.7765252702  71.7335499490 -21.6618891983  19.4690335132  15.1359404176 -73.7530160179  43.5590375474  33.2740457855  -8.0931784795 -39.9709183822 -73.3685716644  74.6665581824  47.3279643342 -44.5065725961 -31.8126048915  -4.1754341347   7.6332146877  -3.3648912920 -14.2235101744 -13.2435315104  72.3224704128 -31.2105340442 -68.0470702542  47.0184569535]);
      PXM=transpose([-0.0120003134   -1.1664196457    4.1530397385   -5.7445045659    2.5902701316    0.7258862806  -46.4317606404  230.9877940490  182.4964553918  -40.6106166252    2.0757079123  -93.1855413667   38.2950387028   94.2575317923    9.2798132122 -380.2771254667  -20.0630096219   22.8324400182   53.2376847335  -20.9514373801  -16.0109579613    2.0165582878  219.1335097025  -86.2136892297   -8.3858784652   -1.2781104875  105.1691056181 -344.8534280980   52.4240278626   39.6526294674]);
      PYm=transpose([4.0411623092     3.9535190774   124.5891817676   -76.9259801505   -86.7355826809     1.3407391280    -6.8391971822  1816.4715564786   751.2401850617   559.9746942482     9.0363178812   -16.5892398645    34.5672651570    -6.0713873958     0.5264778823 -2360.2764269688 -1910.5046325680    -8.2982055205  1214.3992499579    15.7325165165   -12.9666418833    -9.4224851771    81.8388270729   -77.0650114178     1.3566123142    -7.5612094532   139.8405342819  -218.4618272550   -99.0174146424   184.4172020627]);
      PYM=transpose([0.0380243472    0.0658613887  -13.5111317927   13.0900567222   13.1639383905    0.1055632511   -0.9399533190 -573.3268857648 -456.0674353453 -149.2366064701   -3.1759219927   27.4323742558  -14.8932410733   -9.1042120529   -0.0638032940 1017.0431859070  508.7844336416  -54.4047297030 -435.1363541691   40.8572543772   27.1631743937    4.2582421355  -31.5834935503    9.7634205603   45.0627945984    1.5597242582   -7.1062863531  -11.1850871675  -15.1747323477   52.0752889360]);
    case 33
      PXm=transpose([0.1981513061  -2.7514747020   3.8938914121  -2.5478401742  -2.7870179154   2.2798826877  -4.5754816478  58.9996919580 -33.6868807753  13.6147100816  14.5748249749 -82.6649793625  48.8373363446  39.4245991389  -8.1065804413 -14.8208374994 -49.7263854362  68.9073161800  29.3796798494 -40.1569931456 -29.4965853942  -3.4458490773   5.2913974258   0.7302881862 -17.7250530409 -12.8518189342  73.6572238793 -27.9244084624 -70.6081229142  45.4670857196]);
      PXM=transpose([-0.0961727215   -1.2862772214   -0.5101952638   -3.9146770696    4.4577876062    0.8703768827  -44.7859101023  154.1011235551  181.0755931542  -49.0374262874    1.9181723551  -96.8119707011   43.9968285414   91.4390474598    8.9693731996 -314.5786146923   43.5329527291   28.0164971384   17.4657751824  -24.6123038164  -17.6918662006   11.4663088508  173.5425192890  -51.2763425710  -12.3623127004   -1.1197829381   45.3799559665 -237.6733420089   70.9285748194  -24.0975909764]);
      PYm=transpose([4.1272332507     2.9402589806   137.0607533085   -87.6921737865   -89.7787345956     0.3753846514    -9.1049820003  1984.7830822645   806.6561269595   591.1546817938     7.8386254127   -16.5149890247    30.4131620425     4.0843196411    -0.4174661871 -2564.3750075531 -2090.7634906630   -13.1052756930  1346.8104754135    15.9601666935    -4.0449209556    -6.7977191650    72.1553634874   -66.8704779328     0.3645880838    -6.8466124322   130.8126496969  -214.1344282572   -94.2071597850   182.7947200482]);
      PYM=transpose([0.4181946083   -0.3239121491   -9.0442285847   11.5041367031   12.4681381813   -1.6225340693    5.0920065601 -690.2309947626 -553.9754776403 -185.3465784838   -1.9434456438    2.0617178812   -3.9893707435   -6.4343267752   -0.7816178591 1241.7114441950  631.3283339033  -71.4392888751 -545.8795874457   53.7134934915   37.8491367980   -7.2392763059  -27.9506900924   25.8070141409   54.0424897730    0.8563888217   30.7107975796  -44.8846372938  -14.4578111189   39.9793469928]);
    case 34
      PXm=transpose([0.1751974058  -3.4972859865   3.1024827823  -3.0606326353  -1.7780023740   2.3640584179  -4.5576097592  46.8463143881 -43.6705465770   8.1263295746  14.6274095116 -89.2374757192  49.3269245793  45.0941670277  -9.2204634734   2.0615838282 -38.5526604764  73.4413791937  30.6192380489 -43.7502011192 -29.7754238838  -5.1546869594   5.1289422106  -4.9680950009 -19.9947659254 -13.1145348631  84.6289605550 -33.5232724423 -80.9110094482  59.9504068968]);
      PXM=transpose([-0.0288055929   -1.2478581085    5.7515523982   -7.8548653788    2.3530294825    1.0282457407  -48.4628406165  256.2787973273  201.2644243569  -36.9087206555    1.9713541692  -81.8124227538   22.0933737487   95.3784629800    9.8490448017 -423.4185499437  -45.7458859650   27.4399302527   80.3917527336  -22.7816821029  -19.1163973064    3.4808039480  215.8751151249  -85.4515293510   -8.4465679700   -1.2242483678   96.8132413129 -343.6669440519   63.4260714140   37.0814181532]);
      PYm=transpose([3.5900128309     4.1071076845   111.4235059008   -69.1742888032   -77.5230250172     0.6033772258    -7.8029385407  1776.0392364623   711.9187491355   528.3691506196     7.9085594308    -1.3966896761    23.1428245953    -9.6488901820     0.0846633719 -2280.7186598932 -1846.4278533585   -25.1734721712  1171.8965840401    25.8238746508    -1.0371465237    -4.9061754589    50.4832204024   -50.8832764931     4.5306259914    -6.5572493586   112.7727977658  -181.8809381400   -88.3699106796   167.1116728862]);
      PYM=transpose([-0.5695467491   -0.9088954312  -15.5542912440   15.5704723808   19.8016799701    0.4915356742    2.5954945055 -615.4299187164 -495.5344214223 -181.0435696051   -4.3818701518   24.8652012106  -16.7071831510   -8.2631699929   -0.7746115492 1109.8602020991  555.0552217636  -59.1689509927 -481.7406739922   43.2450099390   31.4858262810   -0.8545149324  -35.3413669281   20.6006981653   66.0599542775    2.2200359378    1.1806142612  -26.1820738111   -5.7355266289   42.6210113099]);
    case 35
      PXm=transpose([0.2117074295  -2.9765611503   5.4502095057  -4.0033215364  -3.0492254214   2.2850256025  -5.6608625713  98.5152401813 -30.6540010510  18.5235186097  13.5755100618 -63.4103084592  27.8289800002  39.3071373630  -8.4601679192 -54.5882937105 -93.4030106224  71.9296135447  65.8241980959 -43.9667012466 -28.5658130017  -4.7770414229  13.3532146004 -12.4791587018 -16.8202438485 -12.3857973066  77.8172735718 -54.4868479356 -66.2344534648  71.3230193745]);
      PXM=transpose([-0.1350382620   -1.9360443512   -3.4389468036   -2.4273867365    5.4696085473    1.1060904293  -48.2393960220  161.9971199828  188.8673890371  -53.8704795440    1.5507414232 -103.7166724693   39.3048027974  100.5602350438    9.8901420894 -335.2570555112   58.8559293013   28.2211281659   16.9339833940  -23.0017091345  -19.3910242404   13.2956996803  187.7751720825  -56.6076958870  -17.9907686996   -0.8937057392   43.0036418128 -251.5171028030   82.5641531895  -24.5532965904]);
      PYm=transpose([3.4855615959     2.7330445513   101.2321164792   -67.0559838946   -70.5100054163    -0.7192992800    -5.2092241611  1701.0198143928   715.9735075679   478.4982736443     7.1855773731   -29.5604177727    38.1535075020     1.6480266964     0.2217498753 -2246.6302042085 -1691.8139508211   -38.2725379478  1120.3190179384    30.7742412830    10.8860389902    -7.7959356481    60.0656074256   -48.6963824914   -29.5798839421    -6.0385162833   119.8176561898  -170.0609760116   -90.3560771968   151.0174641672]);
      PYM=transpose([-0.7692490224   -0.7031052200   -8.4947332547   11.8636726083   12.2207364957    3.3746215105    1.9463444241 -431.6249846163 -382.2637993245 -120.9030004188   -7.8196435804   34.9272165662  -30.1914597506   -9.6785329950   -0.3167860208  828.0642237549  357.3376965931  -48.0803271818 -333.4116257463   34.0564050175   26.1035688415   -1.9370660278  -28.1131325187   12.8154070117   58.6826740566    3.7018175957    3.8294433533  -35.3080842166   -0.4798921615   44.9553974116]);
    otherwise
      disp(['pb : donnees']),
  end;

function [xMin,xMax,yMin,yMax]=DMX2BB(pts)
  tailleIm=[120;160];
  n=size(pts,2);
  %Initialization :
  Om=tailleIm(1)*tailleIm(2);
  [xMin,xMax,yMin,yMax]=MAX2BB(pts);
  indE=[];
  indF=find(((xMin==pts(1,:))|(pts(1,:)==xMax))|((yMin==pts(2,:))|(pts(2,:)==yMax)));
  indI=setdiff(1:n,indF);
  xMinOpt=xMin; xMaxOpt=xMax; yMinOpt=yMin; yMaxOpt=yMax;
  while(length(indI)>2)

    %Objective Function J
    S=(xMax-xMin+1)*(yMax-yMin+1);
    nI=length(indI);
    if (Om-S~=0)&(S~=0)&(nI~=0)&(nI~=n)
      J1=(n-length(indI))/(Om-S);
      J2=length(indI)/S;
      J=J1/J2;
    else
      J=0;
    end;

    %Is J Optimal ?
    if (length(indE)==0) JOpt=J;
    else 
      if (JOpt>J) 
        JOpt=J;
        xMinOpt=xMin; xMaxOpt=xMax; yMinOpt=yMin; yMaxOpt=yMax;
      end;
    end;

    %Processing the Boundary : 
    Sk=zeros(1,length(indF));
    for k=1:length(indF)
      [xMinF,xMaxF,yMinF,yMaxF]=MAX2BB(pts(:,[indI setdiff(indF,indF(k))]));
      Sk(k)=(yMaxF-yMinF+1)*(xMaxF-xMinF+1);
    end;
    [a,kMeil]=min(Sk);  

  
    %Moving the Boundary of the BB
    indE=[indE indF(kMeil)];
    indNE=setdiff(1:n,indE);
    [xMin,xMax,yMin,yMax]=MAX2BB(pts(:,indNE));
    indF=indNE(find(((xMin==pts(1,indNE))|(pts(1,indNE)==xMax))|((yMin==pts(2,indNE))|(pts(2,indNE)==yMax))));
    indI=setdiff(1:n,[indE indF]);

  end;
  
  xMin=xMinOpt; xMax=xMaxOpt; yMin=yMinOpt; yMax=yMaxOpt;



function [xMin,xMax,yMin,yMax,sigma,J]=BHM2BB(pts)
  tailleIm=[120;160];
  valMoy=@(n)0.0109+0.2977./n-0.3481./n.^2+0.2887./n.^3;
  [xMin,xMax,yMin,yMax]=BTC2BB(pts);
  p=0.7;
  [HN1,HN2,HN3,HN4]=calculHN(pts);
  ptsAd=[pts(1,:)/tailleIm(2);pts(2,:)/tailleIm(1)];
  X0=x2X(xMin,xMax,yMin,yMax,p);
  options=optimset('fminunc');
options=optimset(options,'Display','off','GradObj','on','MaxFunEvals',10);
  [X,fval,exitflag,output,grad,hessian] = fminunc(@(X) fun(X,ptsAd,HN1,HN2,HN3,HN4,tailleIm),X0,options);
  %output,
  [xMin,xMax,yMin,yMax,p]=X2x(X);
  n=size(pts,2);
  H=Hessi(X.',ptsAd,HN1,HN2,HN3,HN4,tailleIm);
  sigma=inv(H(1:4,1:4))/sqrt(n)*valMoy(n);
  J=fun(X,ptsAd,HN1,HN2,HN3,HN4,tailleIm);

function [score,JV]=fun(X,ptsAd,HN1,HN2,HN3,HN4,tailleIm)
  [xMin,xMax,yMin,yMax,p]=X2x(X);
  [ax,bx,ay,by]=conversion(xMin,xMax,yMin,yMax,tailleIm);
  score=fun211C(p,ax,bx,ay,by,ptsAd,HN1,HN2,HN3,HN4);
  JV=gradJ2(X,ptsAd,HN1,HN2,HN3,HN4,tailleIm);

function J=fun211C(p,ax,bx,ay,by,ptsAd,HN1,HN2,HN3,HN4)

  F1=@(x,a,b)(x-a).*(x<=b).*(a<=x)+(b-a)*(x>b);
  F=@(x,y,ax,bx,ay,by,p)(1-p)*x.*y+p/(bx-ax)/(by-ay)*F1(x,ax,bx).*F1(y,ay,by);
  HT1=F(ptsAd(1,:),ptsAd(2,:),ax,bx,ay,by,p);
  HT2=F(1-ptsAd(1,:),ptsAd(2,:),1-bx,1-ax,ay,by,p);
  HT3=F(1-ptsAd(1,:),1-ptsAd(2,:),1-bx,1-ax,1-by,1-ay,p);
  HT4=F(ptsAd(1,:),1-ptsAd(2,:),ax,bx,1-by,1-ay,p);
  J=1/4*mean((HT1-HN1).^2+(HT2-HN2).^2+(HT3-HN3).^2+(HT4-HN4).^2);


function [p,ax,bx,ay,by]=has();
  p=0.5+0.5*rand(1);
  ax=rand(1);
  ay=rand(1);
  bx=rand(1);
  by=rand(1);

function X=x2X(xMin,xMax,yMin,yMax,p)
  X=[xMin xMax yMin yMax p];

function [xMin,xMax,yMin,yMax,p]=X2x(X)
  xMin=X(1); xMax=X(2); yMin=X(3); yMax=X(4); p=X(5);
 

function [pNv,axNv,bxNv,ayNv,byNv,distance]=dep(p,ax,bx,ay,by,ptsAd,HN1,HN2,HN3,HN4,lambda)
  J=Jac(ptsAd,p,ax,ay,bx,by);
  F=@(x,p,a,b)(1-p)*x+p/(b-a).*(x<=b).*(a<=x).*(x-a)+p*(x>b);
  n=size(ptsAd,2);
  HNVect=reshape([HN1;HN2;HN3;HN4],4*n,1);
  FMat=[F(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by);(ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F(ptsAd(2,:),p,ay,by);(ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*(ones(1,n)-F(ptsAd(2,:),p,ay,by));F(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))];
  FVect=reshape(FMat,4*n,1);
  %dp=(J.'*J+lambda*diag(diag(J.'*J)))\(J.'*(HNVect-FVect));
  dp=(J.'*J+lambda*eye(5))\(J.'*(HNVect-FVect));
  pNv=p+dp(1);
  axNv=ax+dp(2);
  bxNv=bx+dp(3);
  ayNv=ay+dp(4);
  byNv=by+dp(5);
  distance=sqrt(sum(dp(2:5).^2)/4);

function test=estFaux(p,ax,bx,ay,by)
  test1=((0<ax)&&(ax<bx))&&(bx<1);
  test2=((0<ay)&&(ay<by))&&(by<1);
  test=~(test1&&test2);

function [ax,bx,ay,by]=conversion(xMin,xMax,yMin,yMax,tailleIm)
  ax=xMin/tailleIm(2); bx=xMax/tailleIm(2);
  ay=yMin/tailleIm(1); by=yMax/tailleIm(1);

function [p,ax,bx,ay,by,score]=actualiser(pNv,axNv,bxNv,ayNv,byNv,scoreNv);
  score=scoreNv;
  p=pNv;
  ax=axNv;
  ay=ayNv;
  bx=bxNv;
  by=byNv;


function J=Jac(ptsAd,p,ax,ay,bx,by)
%calcul du jacobien
n=size(ptsAd,2);
%init
F=@(x,p,a,b)(1-p)*x+p./(b-a).*(x<=b).*(a<=x).*(x-a)+p*(x>b);
F_p=@(x,p,a,b)-x+(x-a)/(b-a).*(x<=b).*(a<=x)+(x>b); %ok
F_a=@(x,p,a,b)p*(x-b)./(b-a).^2.*(a<=x); %ok
F_b=@(x,p,a,b)-p*(x-a)./(b-a).^2.*(x<=b); %ok
J=zeros(4*n,5);
J(1:4:4*n,1)=(F_p(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by)).'+(F(ptsAd(1,:),p,ax,bx).*F_p(ptsAd(2,:),p,ay,by)).';
J(1:4:4*n,2)=(F_a(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by)).';
J(1:4:4*n,3)=(F_b(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by)).';
J(1:4:4*n,4)=(F(ptsAd(1,:),p,ax,bx).*F_a(ptsAd(2,:),p,ay,by)).';
J(1:4:4*n,5)=(F(ptsAd(1,:),p,ax,bx).*F_b(ptsAd(2,:),p,ay,by)).';
J(2:4:4*n,1)=-(F_p(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by)).'+((ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F_p(ptsAd(2,:),p,ay,by)).';
J(2:4:4*n,2)=-(F_a(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by)).';
J(2:4:4*n,3)=-(F_b(ptsAd(1,:),p,ax,bx).*F(ptsAd(2,:),p,ay,by)).';
J(2:4:4*n,4)=((ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F_a(ptsAd(2,:),p,ay,by)).';
J(2:4:4*n,5)=((ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F_b(ptsAd(2,:),p,ay,by)).';
J(3:4:4*n,1)=-(F_p(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))).'-((ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F_p(ptsAd(2,:),p,ay,by)).';
J(3:4:4*n,2)=-(F_a(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))).';
J(3:4:4*n,3)=-(F_b(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))).';
J(3:4:4*n,4)=-((ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F_a(ptsAd(2,:),p,ay,by)).';
J(3:4:4*n,5)=-((ones(1,n)-F(ptsAd(1,:),p,ax,bx)).*F_b(ptsAd(2,:),p,ay,by)).';
J(4:4:4*n,1)=(F_p(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))).'-(F(ptsAd(1,:),p,ax,bx).*F_p(ptsAd(2,:),p,ay,by)).';
J(4:4:4*n,2)=(F_a(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))).';
J(4:4:4*n,3)=(F_b(ptsAd(1,:),p,ax,bx).*(ones(1,n)-F(ptsAd(2,:),p,ay,by))).';
J(4:4:4*n,4)=-(F(ptsAd(1,:),p,ax,bx).*F_a(ptsAd(2,:),p,ay,by)).';
J(4:4:4*n,5)=-(F(ptsAd(1,:),p,ax,bx).*F_b(ptsAd(2,:),p,ay,by)).';

function [HN1,HN2,HN3,HN4]=calculHN(pts)
  n=size(pts,2);
  for k=1:n
    HN1(k)=1/n*sum((pts(1,1:n)<=pts(1,k))&(pts(2,1:n)<=pts(2,k)));
    HN2(k)=1/n*sum((pts(1,1:n)>pts(1,k))&(pts(2,1:n)<=pts(2,k)));
    HN3(k)=1/n*sum((pts(1,1:n)>pts(1,k))&(pts(2,1:n)>pts(2,k)));
    HN4(k)=1/n*sum((pts(1,1:n)<=pts(1,k))&(pts(2,1:n)>pts(2,k)));
  end;



function [xMin,xMax,yMin,yMax]=conversion2(ax,bx,ay,by,tailleIm)
  xMin=max(1,ax*tailleIm(2));
  xMax=max(1,bx*tailleIm(2));
  yMin=max(1,ay*tailleIm(1));
  yMax=max(1,by*tailleIm(1));



function H=Hessi(x0,ptsAd,HN1,HN2,HN3,HN4,tailleIm)

% initializations

k=length(x0);
hessian=zeros(k,k);
grdd=zeros(k,1);
eps=6.0554544523933429e-6;

% Computaion of stepsize (dh)

ax0=abs(x0);
for i=1:k
    if x0(i,1)~=0
       dax0(i,1)=x0(i,1)/ax0(i,1);
    else
       dax0(i,1)=1;
    end
end

dh=eps*max([ax0, (1e-2)*ones(length(x0),1)]')'.*dax0;
xdh=x0+dh;
dh=xdh-x0;  % This increases precision slightly
ee=[]; I=eye(k);
for i=1:k
    ee(:,i)=I(:,i).*dh;
end


% Computation of f0=f(x0)

%f0=feval(f,x0,varargin{:});
f0=fun(x0,ptsAd,HN1,HN2,HN3,HN4,tailleIm);
% Compute forward step

for i=1:k
%    grdd(i,1)=feval(f,x0+ee(:,i),varargin{:});
     grdd(i,1)=fun(x0+ee(:,i),ptsAd,HN1,HN2,HN3,HN4,tailleIm);
end

% Compute 'double' forward step

for i=1:k
    for j=i:k
%        hessian(i,j)=feval(f,x0+(ee(:,i)+ee(:,j)),varargin{:});
        hessian(i,j)=fun(x0+(ee(:,i)+ee(:,j)),ptsAd,HN1,HN2,HN3,HN4,tailleIm);
        if i~=j
            hessian(j,i)=hessian(i,j);
        end
    end
end

l=ones(1,k); 
grdd=kron(l,grdd); 

H=(((hessian - grdd) - grdd') + f0) ./ (kron(dh,dh'));

function JV=gradJ(X,ptsAd,HN1,HN2,HN3,HN4,tailleIm)
  n=size(ptsAd,2);
  FS=@(x,a,b)(x-a).*(x<=b).*(a<=x)/(b-a)+(x>b);
  F1=@(x,y,ax,bx,ay,by,p)(1-p)*x.*y+p*FS(x,ax,bx).*FS(y,ay,by);
  F2=@(x,y,ax,bx,ay,by,p)(1-p)*(1-x).*y+p*(1-FS(x,ax,bx)).*FS(y,ay,by);
  F3=@(x,y,ax,bx,ay,by,p)(1-p)*(1-x).*(1-y)+p*(1-FS(x,ax,bx)).*(1-FS(y,ay,by));
  F4=@(x,y,ax,bx,ay,by,p)(1-p)*x.*(1-y)+p*FS(x,ax,bx).*(1-FS(y,ay,by));
  F1_p=@(x,y,ax,bx,ay,by,p)-x.*y+FS(x,ax,bx).*FS(y,ay,by);
  F2_p=@(x,y,ax,bx,ay,by,p)-(1-x).*y+(1-FS(x,ax,bx)).*FS(y,ay,by);
  F3_p=@(x,y,ax,bx,ay,by,p)-(1-x).*(1-y)+(1-FS(x,ax,bx)).*(1-FS(y,ay,by));
  F4_p=@(x,y,ax,bx,ay,by,p)-x.*(1-y)+FS(x,ax,bx).*(1-FS(y,ay,by));
  J_5=@(X)1/n/2*sum(F1_p(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F1(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN1)+...
  F2_p(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F2(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN2)+...
  F3_p(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F3(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN3)+...
  F4_p(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F4(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN4));
  FS_a=@(x,a,b)(x-b).*(a<=x).*(x<=b)/(b-a)^2;
  F1_a=@(x,y,ax,bx,ay,by,p)p*FS_a(x,ax,bx).*FS(y,ay,by);
  F2_a=@(x,y,ax,bx,ay,by,p)-p*FS_a(x,ax,bx).*FS(y,ay,by);
  F3_a=@(x,y,ax,bx,ay,by,p)-p*FS_a(x,ax,bx).*(1-FS(y,ay,by));
  F4_a=@(x,y,ax,bx,ay,by,p)p*FS_a(x,ax,bx).*(1-FS(y,ay,by));
  J_1=@(X)1/n/2/tailleIm(2)*sum(F1_a(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F1(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN1)+...
  F2_a(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F2(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN2)+...
  F3_a(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F3(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN3)+...
  F4_a(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F4(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN4));
  FS_b=@(x,a,b)-(x-a).*(a<=x).*(x<=b)/(b-a)^2;
  F1_bx=@(x,y,ax,bx,ay,by,p)p*FS_b(x,ax,bx).*FS(y,ay,by);
  F2_bx=@(x,y,ax,bx,ay,by,p)-p*FS_b(x,ax,bx).*FS(y,ay,by);
  F3_bx=@(x,y,ax,bx,ay,by,p)-p*FS_b(x,ax,bx).*(1-FS(y,ay,by));
  F4_bx=@(x,y,ax,bx,ay,by,p)p*FS_b(x,ax,bx).*(1-FS(y,ay,by));
  J_2=@(X)1/n/2/tailleIm(2)*sum(F1_bx(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F1(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN1)+...
  F2_bx(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F2(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN2)+...
  F3_bx(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F3(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN3)+...
  F4_bx(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F4(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN4));
  F1_ay=@(x,y,ax,bx,ay,by,p)p*FS(x,ax,bx).*FS_a(y,ay,by);
  F2_ay=@(x,y,ax,bx,ay,by,p)p*(1-FS(x,ax,bx)).*FS_a(y,ay,by);
  F3_ay=@(x,y,ax,bx,ay,by,p)-p*(1-FS(x,ax,bx)).*FS_a(y,ay,by);
  F4_ay=@(x,y,ax,bx,ay,by,p)-p*FS(x,ax,bx).*FS_a(y,ay,by);
  J_3=@(X)1/n/2/tailleIm(1)*sum(F1_ay(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F1(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN1)+...
  F2_ay(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F2(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN2)+...
  F3_ay(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F3(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN3)+...
  F4_ay(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F4(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN4));
  F1_by=@(x,y,ax,bx,ay,by,p)p*FS(x,ax,bx).*FS_b(y,ay,by);
  F2_by=@(x,y,ax,bx,ay,by,p)p*(1-FS(x,ax,bx)).*FS_b(y,ay,by);
  F3_by=@(x,y,ax,bx,ay,by,p)-p*(1-FS(x,ax,bx)).*FS_b(y,ay,by);
  F4_by=@(x,y,ax,bx,ay,by,p)-p*FS(x,ax,bx).*FS_b(y,ay,by);
  J_4=@(X)1/n/2/tailleIm(1)*sum(F1_by(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F1(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN1)+...
  F2_by(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F2(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN2)+...
  F3_by(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F3(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN3)+...
  F4_by(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5)).*(F4(ptsAd(1,:),ptsAd(2,:),X(1)/tailleIm(2),X(2)/tailleIm(2),X(3)/tailleIm(1),X(4)/tailleIm(1),X(5))-HN4));
  JV=[J_1(X); J_2(X); J_3(X); J_4(X); J_5(X)];

function JV=gradJ2(X,ptsAd,HN1,HN2,HN3,HN4,tailleIm)
  n=size(ptsAd,2);
  FS=@(x,a,b)(x-a).*(x<=b).*(a<=x)/(b-a)+(x>b);
  ax=X(1)/tailleIm(2); bx=X(2)/tailleIm(2); ay=X(3)/tailleIm(1); by=X(4)/tailleIm(1); p=X(5);
  F1=@(x,y)(1-p)*x.*y+p*FS(x,ax,bx).*FS(y,ay,by);
  F1V=F1(ptsAd(1,:),ptsAd(2,:));
  F2=@(x,y)(1-p)*(1-x).*y+p*(1-FS(x,ax,bx)).*FS(y,ay,by);
  F2V=F2(ptsAd(1,:),ptsAd(2,:));
  F3=@(x,y)(1-p)*(1-x).*(1-y)+p*(1-FS(x,ax,bx)).*(1-FS(y,ay,by));
  F3V=F3(ptsAd(1,:),ptsAd(2,:));
  F4=@(x,y)(1-p)*x.*(1-y)+p*FS(x,ax,bx).*(1-FS(y,ay,by));
  F4V=F4(ptsAd(1,:),ptsAd(2,:));
  F1_p=@(x,y)-x.*y+FS(x,ax,bx).*FS(y,ay,by);
  F2_p=@(x,y)-(1-x).*y+(1-FS(x,ax,bx)).*FS(y,ay,by);
  F3_p=@(x,y)-(1-x).*(1-y)+(1-FS(x,ax,bx)).*(1-FS(y,ay,by));
  F4_p=@(x,y)-x.*(1-y)+FS(x,ax,bx).*(1-FS(y,ay,by));
  J_5=1/n/2*sum(F1_p(ptsAd(1,:),ptsAd(2,:)).*(F1V-HN1)+...
  F2_p(ptsAd(1,:),ptsAd(2,:)).*(F2V-HN2)+...
  F3_p(ptsAd(1,:),ptsAd(2,:)).*(F3V-HN3)+...
  F4_p(ptsAd(1,:),ptsAd(2,:)).*(F4V-HN4));
  FS_a=@(x,a,b)(x-b).*(a<=x).*(x<=b)/(b-a)^2;
  F1_a=@(x,y)p*FS_a(x,ax,bx).*FS(y,ay,by);
  F2_a=@(x,y)-p*FS_a(x,ax,bx).*FS(y,ay,by);
  F3_a=@(x,y)-p*FS_a(x,ax,bx).*(1-FS(y,ay,by));
  F4_a=@(x,y)p*FS_a(x,ax,bx).*(1-FS(y,ay,by));
  J_1=1/n/2/tailleIm(2)*sum(F1_a(ptsAd(1,:),ptsAd(2,:)).*(F1V-HN1)+...
  F2_a(ptsAd(1,:),ptsAd(2,:)).*(F2V-HN2)+...
  F3_a(ptsAd(1,:),ptsAd(2,:)).*(F3V-HN3)+...
  F4_a(ptsAd(1,:),ptsAd(2,:)).*(F4V-HN4));
  FS_b=@(x,a,b)-(x-a).*(a<=x).*(x<=b)/(b-a)^2;
  F1_bx=@(x,y)p*FS_b(x,ax,bx).*FS(y,ay,by);
  F2_bx=@(x,y)-p*FS_b(x,ax,bx).*FS(y,ay,by);
  F3_bx=@(x,y)-p*FS_b(x,ax,bx).*(1-FS(y,ay,by));
  F4_bx=@(x,y)p*FS_b(x,ax,bx).*(1-FS(y,ay,by));
  J_2=1/n/2/tailleIm(2)*sum(F1_bx(ptsAd(1,:),ptsAd(2,:)).*(F1V-HN1)+...
  F2_bx(ptsAd(1,:),ptsAd(2,:)).*(F2V-HN2)+...
  F3_bx(ptsAd(1,:),ptsAd(2,:)).*(F3V-HN3)+...
  F4_bx(ptsAd(1,:),ptsAd(2,:)).*(F4V-HN4));
  F1_ay=@(x,y)p*FS(x,ax,bx).*FS_a(y,ay,by);
  F2_ay=@(x,y)p*(1-FS(x,ax,bx)).*FS_a(y,ay,by);
  F3_ay=@(x,y)-p*(1-FS(x,ax,bx)).*FS_a(y,ay,by);
  F4_ay=@(x,y)-p*FS(x,ax,bx).*FS_a(y,ay,by);
  J_3=1/n/2/tailleIm(1)*sum(F1_ay(ptsAd(1,:),ptsAd(2,:)).*(F1V-HN1)+...
  F2_ay(ptsAd(1,:),ptsAd(2,:)).*(F2V-HN2)+...
  F3_ay(ptsAd(1,:),ptsAd(2,:)).*(F3V-HN3)+...
  F4_ay(ptsAd(1,:),ptsAd(2,:)).*(F4V-HN4));
  F1_by=@(x,y)p*FS(x,ax,bx).*FS_b(y,ay,by);
  F2_by=@(x,y)p*(1-FS(x,ax,bx)).*FS_b(y,ay,by);
  F3_by=@(x,y)-p*(1-FS(x,ax,bx)).*FS_b(y,ay,by);
  F4_by=@(x,y)-p*FS(x,ax,bx).*FS_b(y,ay,by);
  J_4=1/n/2/tailleIm(1)*sum(F1_by(ptsAd(1,:),ptsAd(2,:)).*(F1V-HN1)+...
  F2_by(ptsAd(1,:),ptsAd(2,:)).*(F2V-HN2)+...
  F3_by(ptsAd(1,:),ptsAd(2,:)).*(F3V-HN3)+...
  F4_by(ptsAd(1,:),ptsAd(2,:)).*(F4V-HN4));
  JV=[J_1; J_2; J_3; J_4; J_5];
