How to find major and minor axis of a polygon for polygon subdivision

I want to divide a polygon into smaller polygons. For this purpose my procedure was to first convert the polygon into minimum bounding geometry and divide it into smaller polygons using major and minor axis and finally clip it using the main polygon. I have written a code but it divide the polygon wrongly, this is because for finding the major axis of the minimum bounding geometry, I have to use its envelope, and because the orientation of envelope is different from minimum bounding geometry, the line axis for dividing the polygon is determined wrongly. So, I searched the net and found the same problem that is solved in VBA. I tried to convert the code to C#. The link is as follows: https://stackoverflow.com/questions/15291223/convex-quadrilateral-polygons-subdivision-in-equal-parts-using-python-ogr/16355198#16355198

And my code is as following:

const int xSplit = 5;
const int ySplit = 5;
Double[] dCoords=new double;
int lID;

/////////////////// select features
IPolygon pPolygon = new PolygonClass();
IFeatureCursor polygonFeatCursor = MBRFeatureClass.Search(null, false);
IFeature polygonFeature = polygonFeatCursor.NextFeature();
while (polygonFeature!=null)
{
pPolygon = polygonFeature.Shape as IPolygon;
polygonFeature = polygonFeatCursor.NextFeature();
}

private void GridQuadriLateral(IPolygon pPolygon, IFeatureClass buildingFeatureClass)
{
ISegmentCollection pSegmentCollection = pPolygon as ISegmentCollection;
const int xSplit = 5;
const int ySplit = 5;
Double[ , , ] dCoords =new double[xSplit,ySplit, 2];
int lIdx=0;
double[] cx = new double;
double[] cy = new double;
double[] dx = new double;
double[] dy = new double;
double dx2;
double dy2;
double x1;
double x2;
double x3;
double y1;
double y2;
double y3;
int l;
int xs;
int ys;

// Get the corner coords of the quad
for (l = 0; l <= 3; l++)
{
lIdx = GetIndexOfNextCornerSegment(lIdx, pPolygon);
cx[l] = pSegmentCollection.Segment[l].FromPoint.X;
cy[l] = pSegmentCollection.Segment[l].FromPoint.Y;
}

dx = ((cx - cx)/xSplit);
dx = ((cx - cx) / ySplit);
dx = ((cx - cx) / xSplit);
dx = ((cx - cx) / ySplit);

dy = ((cy - cy) / xSplit);
dy = ((cy - cy) / ySplit);
dy = ((cy - cy) / xSplit);
dy = ((cy - cy) / ySplit);

for (ys= 0; ys<ySplit; ys++)
{
x1 = cx + dx * ys;
y1 = cy + dy * ys;
x2 = cx + dx * ys;
y2 = cy + dy * ys;
dx2 = ((x2 - x1) / xSplit);
dy2 = ((y2 - y1) / ySplit);

for (xs = 0; xs <xSplit; xs++)
{
x3 = x1 + dx2 * xs;
y3 = y1 + dy2 * ys;
dCoords[xs, ys,0] = x3;
dCoords[xs, ys,1] = y3;

}
}

//build grid
IFeatureBuffer pFtrBfr = null;
pFtrBfr = buildingFeatureClass.CreateFeatureBuffer();
IPointCollection pPointCollection ;
IPoint pPoint;
IFeatureCursor pFeatureCursor = buildingFeatureClass.Insert(true);
for (int i = 0; i < ySplit-1; i++)
{
for (int j = 0; j < xSplit-1; j++)
{
pPointCollection = new PolygonClass();
pPoint = new PointClass();
pPoint.PutCoords(dCoords[j,i,0], dCoords[j,i,1]);
pPoint.PutCoords(dCoords[j, i+1, 0], dCoords[j, i+1, 1]);
pPoint.PutCoords(dCoords[j+1, i+1, 0], dCoords[j+1, i+1, 1]);
pPoint.PutCoords(dCoords[j+1, i, 0], dCoords[j+1, i, 1]);
pPoint.PutCoords(dCoords[j, i, 0], dCoords[j, i, 1]);

pFtrBfr.Shape = pPointCollection as IGeometry;
pFeatureCursor.InsertFeature(pFtrBfr);
pFeatureCursor.Flush();

/*
IPolygon ppPolygon = new PolygonClass();
ppPolygon = pPointCollection as IPolygon;
ppPolygon.Close();

IFeature pFeat = buildingFeatureClass.CreateFeature();
pFeat.Shape = ppPolygon;
pFeat.Store();
*/

}
}

m_hookHelper.ActiveView.Refresh();
//throw new NotImplementedException();
}

private int GetIndexOfNextCornerSegment(int lStartIdx, IPolygon pPolygon)
{
int index = 0;

ISegmentCollection pSegmentCollection = pPolygon as ISegmentCollection;
ILine pLine1 = new LineClass();
ILine pLine2 = new LineClass();
int l;
int lNxtIdx;
double dAng;

for ( l = 0; l < pSegmentCollection.SegmentCount-1; l++)
{
lNxtIdx = lStartIdx + 1;
if (lNxtIdx==pSegmentCollection.SegmentCount)
{
lNxtIdx = 0;
}
pLine1 = pSegmentCollection.Segment[lStartIdx] as ILine;
lNxtIdx = lNxtIdx + 1;
if (lNxtIdx==pSegmentCollection.SegmentCount)
{
lNxtIdx = 0;
}
pLine2 = pSegmentCollection.Segment[lNxtIdx] as ILine;
dAng = Math.Abs(pLine1.Angle-pLine2.Angle)*(180/(Math.PI));
if (dAng>=180)
{
dAng = 360 - dAng;
}
if (dAng>20)
{
index= lNxtIdx;
break;
}
}

return index;
//throw new NotImplementedException();
}

But now my problem is that when I run the code, the output is not partitioned completely. I am completely confused about it. The image of the output is as follows: • A rectangle is the trivial case. All you need is opposite midpoints and atan2. – Vince May 30 '18 at 12:24
• @Vince Thank you very much for your comment. Can you guide me how to use it for the definition of the line? – parastoo May 30 '18 at 12:59
• No, I don't think I can. Certainly not from a cell phone. Plot your polygon on a piece of paper, tic off midpoints on the four segments, and connect the opposites. Then do it again in code. You can toss all the MBR calculation code, because it's not of any use. – Vince May 30 '18 at 14:11
• For a more general case you'll need to compute the convex hull and then use an algorithm, such as rotating calipers, to compute the minimum bounding box. – danielm May 30 '18 at 21:49
• – Kirk Kuykendall May 31 '18 at 3:42

Finally, I could partition a minimum bounding rectangle into smaller units using the following code:

IPolygon pPolygon = new PolygonClass();
IFeatureCursor polygonFeatCursor = MBRFeatureClass.Search(null, false);
IFeature polygonFeature = polygonFeatCursor.NextFeature();
while (polygonFeature!=null)
{
pPolygon = polygonFeature.Shape as IPolygon;
int ID = Convert.ToInt32(polygonFeature.get_Value(polygonFeature.Fields.FindField("OBJECTID")));

polygonFeature.Delete();
polygonFeature = polygonFeatCursor.NextFeature();
}

MessageBox.Show("Finished!!!!!!!!!!!");

private void GridQuadriLateral(IPolygon pPolygon, IFeature polygonFeature, IFeatureClass MBRFeatureClass)
{
double length=10;
double width=10;
double[] cx = new double;
double[] cy = new double;

ISegmentCollection pSegmentCollection = pPolygon as ISegmentCollection;
// Get the corner coords of the quad
for (int l = 0; l <= 3; l++)
{
cx[l] = pSegmentCollection.Segment[l].FromPoint.X;
cy[l] = pSegmentCollection.Segment[l].FromPoint.Y;
}

IPoint firstPoint = new PointClass();
IPoint zeroPoint = new PointClass();
IPoint mainZeroPoint = new PointClass();
double mainAngle = Convert.ToDouble(polygonFeature.get_Value(polygonFeature.Fields.FindField("MBG_Orient")));
#region main angle >=90
if (mainAngle>=90)
{
//point2 is the point with lowest Y value
firstPoint.X = cx;
firstPoint.Y = cy;
double x = cx - cx;
double y = cy - cy;
double s = (x * x) + (y * y);
double hypot = Math.Sqrt(s); //this is the length of the baseline
double myAngle = (Math.Atan(y / x)) * (180 / Math.PI);

double p = cx - cx;
double q = cy - cy;
double t = (p * p) + (q * q);
double hypot2 = Math.Sqrt(t);  //this is the height or vertical dist
double myAngle2 = (Math.Atan(q / p)) * (180 / Math.PI);
///////////////////////////////////////
if (hypot<=length*1.5 && hypot2<=length*1.5)
{
return;
}
else
{
// find how many small units to create

int numDiv2 = Convert.ToInt32(hypot2/width);
int numDiv1 = Convert.ToInt32(hypot / length);
double dist2 = Convert.ToDouble(hypot2/numDiv2);
double dist1 = Convert.ToDouble(hypot / numDiv1);

//double p1 = dist1;
Double[, ,] myArray = new double[numDiv1 * numDiv2, 4, 2];
IPoint basePoint1 = new PointClass();
IPoint horPoint1 = new PointClass();
IPoint ePoint1 = new PointClass();
List<ILine> lineList1 = new List<ILine>();
IPointCollection pPointCollection;
IFeatureBuffer pFtrBfr = null;
pFtrBfr = MBRFeatureClass.CreateFeatureBuffer();
IFeatureCursor pFeatureCursor = MBRFeatureClass.Insert(true);
zeroPoint = firstPoint;
mainZeroPoint = firstPoint;
for (int i = 0; i < numDiv2; i++)
{
zeroPoint = GetPoint(mainZeroPoint, myAngle2, dist2 * i);
for (int j = 0; j < numDiv1; j++)
{
firstPoint = GetPoint(zeroPoint, myAngle, dist1 * (j));
pPointCollection = new PolygonClass();
myArray[i, 0, 0] = firstPoint.X;
myArray[i, 0, 1] = firstPoint.Y;
basePoint1 = GetPoint(firstPoint, myAngle2, dist2);
myArray[i, 1, 0] = basePoint1.X;
myArray[i, 1, 1] = basePoint1.Y;
horPoint1 = GetPoint(basePoint1, myAngle, dist1);
myArray[i, 2, 0] = horPoint1.X;
myArray[i, 2, 1] = horPoint1.Y;
ePoint1 = GetPoint(firstPoint, myAngle, dist1);
myArray[i, 3, 0] = ePoint1.X;
myArray[i, 3, 1] = ePoint1.Y;

pFtrBfr.Shape = pPointCollection as IGeometry;
pFeatureCursor.InsertFeature(pFtrBfr);
pFeatureCursor.Flush();
}
}
}
}
#endregion

#region main angle < 90
if (mainAngle<90)
{
//point3 is the point with lowest Y value
firstPoint.X = cx;
firstPoint.Y = cy;

double x = cx - cx;
double y = cy - cy;
double s = (x * x) + (y * y);
double hypot = Math.Sqrt(s); //this is the length of the baseline
double myAngle = (Math.Atan(y / x)) * (180 / Math.PI);

double p = cx - cx;
double q = cy - cy;
double t = (p * p) + (q * q);
double hypot2 = Math.Sqrt(t);  //this is the height or vertical dist
double myAngle2 = (Math.Atan(q / p)) * (180 / Math.PI);
///////////////////////////////////////
if (hypot <= length * 1.5 && hypot2 <= length * 1.5)
{
return;
}
else
{
// find how many small units to create

int numDiv2 = Convert.ToInt32(hypot2 / width);
int numDiv1 = Convert.ToInt32(hypot / length);
double dist2 = Convert.ToDouble(hypot2 / numDiv2);
double dist1 = Convert.ToDouble(hypot / numDiv1);

//double p1 = dist1;
Double[, ,] myArray = new double[numDiv1 * numDiv2, 4, 2];
IPoint basePoint1 = new PointClass();
IPoint horPoint1 = new PointClass();
IPoint ePoint1 = new PointClass();
List<ILine> lineList1 = new List<ILine>();
IPointCollection pPointCollection;
IFeatureBuffer pFtrBfr = null;
pFtrBfr = MBRFeatureClass.CreateFeatureBuffer();
IFeatureCursor pFeatureCursor = MBRFeatureClass.Insert(true);
zeroPoint = firstPoint;
mainZeroPoint = firstPoint;
for (int i = 0; i < numDiv2; i++)
{
zeroPoint = GetPoint(mainZeroPoint, myAngle2, dist2 * i);
for (int j = 0; j < numDiv1; j++)
{
firstPoint = GetPoint(zeroPoint, myAngle, dist1 * (j));
pPointCollection = new PolygonClass();
myArray[i, 0, 0] = firstPoint.X;
myArray[i, 0, 1] = firstPoint.Y;
basePoint1 = GetPoint(firstPoint, myAngle2, dist2);
myArray[i, 1, 0] = basePoint1.X;
myArray[i, 1, 1] = basePoint1.Y;
horPoint1 = GetPoint(basePoint1, myAngle, dist1);
myArray[i, 2, 0] = horPoint1.X;
myArray[i, 2, 1] = horPoint1.Y;
ePoint1 = GetPoint(firstPoint, myAngle, dist1);
myArray[i, 3, 0] = ePoint1.X;
myArray[i, 3, 1] = ePoint1.Y;

pFtrBfr.Shape = pPointCollection as IGeometry;
pFeatureCursor.InsertFeature(pFtrBfr);
pFeatureCursor.Flush();

}
}
}

}

#endregion

m_hookHelper.ActiveView.Refresh();
}

private IPoint GetPoint(IPoint basePoint1, double myAngle, double hypot)
{
double distance=hypot;
IPoint endPoint = new PointClass();
if (myAngle<0)
{
myAngle = myAngle + 180;
}
double offsetX = distance * Math.Cos((myAngle) * (Math.PI / 180));
double offsetY = distance * Math.Sin((myAngle) * (Math.PI / 180));
endPoint.X = (basePoint1.X) + (offsetX);
endPoint.Y = (basePoint1.Y) + (offsetY);
return endPoint;
//throw new NotImplementedException();
}