Polygon Subdivision – Finding Major and Minor Axis

arcobjectscpolygon

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[3];
            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;
                GridQuadriLateral(pPolygon, buildingFeatureClass);
                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[4];
            double[] cy = new double[4];
            double[] dx = new double[4];
            double[] dy = new double[4];
            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[0] = ((cx[1] - cx[0])/xSplit);
            dx[1] = ((cx[1] - cx[2]) / ySplit);
            dx[2] = ((cx[2] - cx[3]) / xSplit);
            dx[3] = ((cx[0] - cx[3]) / ySplit);

            dy[0] = ((cy[1] - cy[0]) / xSplit);
            dy[1] = ((cy[1] - cy[2]) / ySplit);
            dy[2] = ((cy[2] - cy[3]) / xSplit);
            dy[3] = ((cy[0] - cy[3]) / ySplit);

            for (ys= 0; ys<ySplit; ys++)
            {
                x1 = cx[3] + dx[3] * ys;
                y1 = cy[3] + dy[3] * ys;
                x2 = cx[2] + dx[1] * ys;
                y2 = cy[2] + dy[1] * 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]);
                    pPointCollection.AddPoint(pPoint);
                    pPoint.PutCoords(dCoords[j, i+1, 0], dCoords[j, i+1, 1]);
                    pPointCollection.AddPoint(pPoint);
                    pPoint.PutCoords(dCoords[j+1, i+1, 0], dCoords[j+1, i+1, 1]);
                    pPointCollection.AddPoint(pPoint);
                    pPoint.PutCoords(dCoords[j+1, i, 0], dCoords[j+1, i, 1]);
                    pPointCollection.AddPoint(pPoint);
                    pPoint.PutCoords(dCoords[j, i, 0], dCoords[j, i, 1]);
                    pPointCollection.AddPoint(pPoint);

                    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:
enter image description here

Best Answer

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")));

                    GridQuadriLateral(pPolygon, polygonFeature, MBRFeatureClass);

                    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[4];
        double[] cy = new double[4];

        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[2];
            firstPoint.Y = cy[2];
            double x = cx[1] - cx[2];
            double y = cy[1] - cy[2];
            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[0] - cx[1];
            double q = cy[0] - cy[1];
            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;
                        pPointCollection.AddPoint(firstPoint);
                        basePoint1 = GetPoint(firstPoint, myAngle2, dist2);
                        myArray[i, 1, 0] = basePoint1.X;
                        myArray[i, 1, 1] = basePoint1.Y;
                        pPointCollection.AddPoint(basePoint1);
                        horPoint1 = GetPoint(basePoint1, myAngle, dist1);
                        myArray[i, 2, 0] = horPoint1.X;
                        myArray[i, 2, 1] = horPoint1.Y;
                        pPointCollection.AddPoint(horPoint1);
                        ePoint1 = GetPoint(firstPoint, myAngle, dist1);
                        myArray[i, 3, 0] = ePoint1.X;
                        myArray[i, 3, 1] = ePoint1.Y;
                        pPointCollection.AddPoint(ePoint1);
                        pPointCollection.AddPoint(firstPoint);

                        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[3];
            firstPoint.Y = cy[3];

            double x = cx[2] - cx[3];
            double y = cy[2] - cy[3];
            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[1] - cx[2];
            double q = cy[1] - cy[2];
            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;
                        pPointCollection.AddPoint(firstPoint);
                        basePoint1 = GetPoint(firstPoint, myAngle2, dist2);
                        myArray[i, 1, 0] = basePoint1.X;
                        myArray[i, 1, 1] = basePoint1.Y;
                        pPointCollection.AddPoint(basePoint1);
                        horPoint1 = GetPoint(basePoint1, myAngle, dist1);
                        myArray[i, 2, 0] = horPoint1.X;
                        myArray[i, 2, 1] = horPoint1.Y;
                        pPointCollection.AddPoint(horPoint1);
                        ePoint1 = GetPoint(firstPoint, myAngle, dist1);
                        myArray[i, 3, 0] = ePoint1.X;
                        myArray[i, 3, 1] = ePoint1.Y;
                        pPointCollection.AddPoint(ePoint1);
                        pPointCollection.AddPoint(firstPoint);

                        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();
            }
Related Question