在python中如何、何时以及什么样的向量化？

2024-10-04 05:26:09 发布

您现在位置：Python中文网/ 问答频道 /正文
9091 3

对，所以这基本上是我前面一个问题的后续。我有一些二进制数据是浮点二进制格式的。使用C时，该过程很快，但使用atof（）会损失一些精度。我试着浏览论坛，还有其他地方，但我的问题没有得到解决。因此，我搬到了python。啊，快乐！这个程序运行得非常好，但与C相比，它的运行速度非常慢。我在python上找到了优化，这让我想到了Cython和Weave，但我有一些疑问。如果您按照我的代码进行操作，我很困惑在哪里应用优化的C代码，因为我是从numpy对象读取的。我的问题是，是否可以在Cython中使用numpy函数读取数据，如果可以，请提供一个小例子。在

C代码使用PolSARpro的头文件，libbmp用于创建.bmp文件

作为一个注释，我把我的代码都贴出来了。天知道，为了让公式发挥作用，我不得不经历很多。其他人也可以这样做：）

C代码（正常工作，但atof（）会失去精度，因此输出lat long有点不正常）

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <polSARpro/bmpfile.c>
#include <polSARpro/graphics.c>
#include <polSARpro/matrix.c>
#include <polSARpro/processing.c>
#include <polSARpro/util.c>
#define METAL_THRESHOLD 5.000000
#define POLARIZATION_FRACTION_THRESHOLD 0.900000
#define PI 3.14159265
#define FOURTHPI PI/4
#define deg2rad PI/180
#define rad2deg 180./PI

/*double PI = 3.14159265;
double FOURTHPI = PI / 4;
double deg2rad = PI / 180;
double rad2deg = 180.0 / PI;*/

FILE *L1,*PF,*SPF;
FILE *txt;
FILE *finalLocations;
long i=0,loop_end;
int lig,col;
float l1,pf,spf;
long pos;
int Nlig,Ncol;

float *bufferout;
float *bufferin_L1,*bufferin_L2;
float valueL1,valuePF,xx;
float sizeGridX, sizeGridY, startX, startY;
float posX,posY;
int ZONE;
char Heading[10];
char setZone[15];

int p[4][2];

int degree, minute, second;

void UTM2LL(int ReferenceEllipsoid, double UTMNorthing, double UTMEasting, char* UTMZone, double *Lat, double *Long)
{
//converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees.
    //Written by Chuck Gantz- chuck.gantz@globalstar.com

    double k0 = 0.9996;
    double a = 6378137;
    double eccSquared =  0.00669438;
    double eccPrimeSquared;
    double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
    double N1, T1, C1, R1, D, M;
    double LongOrigin;
    double mu, phi1, phi1Rad;
    double x, y;
    int ZoneNumber;
    char* ZoneLetter;
    int NorthernHemisphere; //1 for northern hemispher, 0 for southern

    x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
    y = UTMNorthing;

    ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10);
    if((*ZoneLetter - 'N') >= 0)
        NorthernHemisphere = 1;//point is in northern hemisphere
    else
    {
        NorthernHemisphere = 0;//point is in southern hemisphere
        y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
    }

    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone

    eccPrimeSquared = (eccSquared)/(1-eccSquared);

    M = y / k0;
    mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));

    phi1Rad = mu    + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
                + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
                +(151*e1*e1*e1/96)*sin(6*mu);
    phi1 = phi1Rad*rad2deg;

    N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
    T1 = tan(phi1Rad)*tan(phi1Rad);
    C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
    R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
    D = x/(N1*k0);

    *Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
                    +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
    *Lat = *Lat * rad2deg;

    *Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
                    *D*D*D*D*D/120)/cos(phi1Rad);
    *Long = LongOrigin + *Long * rad2deg;
}

void convertToDegree(float decimal)
{
    int negative = decimal < 0;
    decimal = abs(decimal);
    minute = (decimal * 3600/ 60);
    second = fmodf((decimal * 3600),60);
    degree = minute / 60;
    minute = minute % 60;
    if (negative)
    {
        if (degree > 0)
            degree = -degree;
        else if (minute > 0)
            minute = -minute;
        else
            second = -second;
    }
}

void readConfig(int *Row, int *Col)
{
    char tmp[70];
    int i=0;
    FILE *fp = fopen("config.txt","r");
    if(fp == NULL)
    {
        perror("Config.txt");
        exit(1);
    }
    while(!feof(fp))
    {
        fgets(tmp,70,fp);
        if (i==1)
            *Row = atoi(tmp);
        if(i==4)
            *Col = atoi(tmp);
        i++;
    }
    fclose(fp);
}

void readHDR(float *gridX,float *gridY,float *startXPos,float *startYPos)
{
    FILE *fp = fopen("PF.bin.hdr","r");
    int i=0;
    char tmp[255];
    char junk[255];
    memset(tmp,0X00,sizeof(tmp));
    memset(junk,0X00,sizeof(junk));
    if(fp==NULL)
    {
        perror("Please locate or create PF.bin.hdr");
        exit(0);
    }
    while(!feof(fp))
    {
        if(i==13)
            break;
        fgets(tmp,255,fp);
        i++;
    }
    fclose(fp);

    strcpy(junk,strtok(tmp,","));
    strtok(NULL,",");
    strtok(NULL,",");
    strcpy(tmp,strtok(NULL,","));
    //puts(tmp);
    *startXPos = atof(tmp);
    strcpy(tmp,strtok(NULL,","));
    //puts(tmp);
    *startYPos = atof(tmp);
    strcpy(tmp,strtok(NULL,","));
    //puts(tmp);
    *gridX = atof(tmp);
    strcpy(tmp,strtok(NULL,","));
    //puts(tmp);
    *gridY = atof(tmp);
    strcpy(tmp,strtok(NULL,","));
    ZONE = atoi(tmp);
    strcpy(tmp,strtok(NULL,","));
    strcpy(Heading,tmp);
}

int main()
{
    bmpfile_t *bmp;
    double Lat;
    double Long;
    int i;
    rgb_pixel_t pixelMetal = {128, 64, 0, 0};
    rgb_pixel_t pixelOthers = {128, 64, 0, 0};
    readConfig(&Nlig,&Ncol);
    readHDR(&sizeGridX,&sizeGridY,&startX,&startY);
    //startX = startX - (double) 0.012000;
    //startY = startY + (double)0.111000;
    printf("Enter the rectangle's top-left and bottom-right region of interest points as: x y\n");
    for(i=0;i<2;i++)
    {
        printf("Enter point %d::\t",i+1);
        scanf("%d %d",&p[i][0], &p[i][1]);
    }
    printf("Grid Size(X,Y)::( %f,%f ), Start Positions(X,Y)::( %f, %f ), ZONE::%d, Heading:: %s\n\n",sizeGridX,sizeGridY,startX,startY,ZONE,Heading);
    pixelMetal.red = 255;
    pixelMetal.blue = 010;
    pixelMetal.green = 010;
    pixelOthers.red = 8;
    pixelOthers.blue = 8;
    pixelOthers.green = 8;
    L1 = fopen("l1.bin","rb");
    PF =fopen("PF.bin","rb");
    SPF = fopen("SPF_L1.bin","wb");
    //txt = fopen("locations(UTM).txt","w");
    finalLocations = fopen("locationsROI.txt","w");
    if(L1==NULL || PF==NULL || SPF==NULL || finalLocations == NULL)
    {
        perror("Error in opening files!");
        return -1;
    }
    fseek(L1,0,SEEK_END);
    pos = ftell(L1);
    loop_end = pos;
    printf("L1.bin contains::\t%ld elements\n",pos);
    fseek(PF,0,SEEK_END);
    pos = ftell(PF);
    printf("PF.bin contains::\t%ld elements\n",pos);
    fseek(L1,0,SEEK_SET);
    fseek(PF,0,SEEK_SET);
    bmp = bmp_create(Ncol,Nlig,8); //width * height
    bufferin_L1 = vector_float(Ncol);
    bufferin_L2 = vector_float(Ncol);
    bufferout = vector_float(Ncol);
    printf("Resources Allocated. Beginning...\n");
    for (lig = 0; lig < Nlig; lig++) /* rows */
    {
        if (lig%(int)(Nlig/20) == 0)
        {
            printf("%f\r", 100. * lig / (Nlig - 1));
            fflush(stdout);
        }
        fread(&bufferin_L1[0], sizeof(float), Ncol, L1);
        fread(&bufferin_L2[0], sizeof(float), Ncol, PF);
        for (col = 0; col < Ncol; col++) /* columns */
        {
            valueL1 = bufferin_L1[col];
            valuePF = bufferin_L2[col];
            if(valueL1 >= METAL_THRESHOLD && valuePF >= POLARIZATION_FRACTION_THRESHOLD)
            {
                if(col >= p[0][0] && col <= p[1][0] && lig >= p[0][1] && lig <= p[1][1])
                {
                    xx = fabs(valueL1 + valuePF);
                    bmp_set_pixel(bmp,col,lig,pixelMetal);
                    posX = startX + (sizeGridX * col);
                    posY = startY - (sizeGridY * lig);
                    //fprintf(txt,"%f %f %d %s\n",posX,posY,ZONE,Heading);
                    sprintf(setZone,"%d",ZONE);
                    if(strstr(Heading,"Nor")!=NULL)
                        strcat(setZone,"N");
                    else
                        strcat(setZone,"S");
                    UTM2LL(23, posY, posX, setZone, &Lat, &Long); // 23 for WGS-84
                    convertToDegree(Lat);
                    //fprintf(finalLocations,"UTM:: %.2fE %.2fN , Decimal: %f %f , Degree: %d %d %d, ",posX,posY,Lat,Long,degree,minute,second);
                    //fprintf(finalLocations,"%.2fE,%.2fN,%f,%f ,%d,%d,%d,",posX,posY,Lat,Long,degree,minute,second);
                    fprintf(finalLocations,"%.2f,%.2f,%f,%f ,%d,%d,%d,",posX,posY,Lat,Long,degree,minute,second);
                    convertToDegree(Long);
                    fprintf(finalLocations,"%d,%d,%d\n",degree,minute,second);
                }
                else
                {
                    xx = fabs(valueL1) ;
                    bmp_set_pixel(bmp,col,lig,pixelOthers);
                }
            }
            else
            {
                xx = fabs(valueL1) ;
                bmp_set_pixel(bmp,col,lig,pixelOthers);
            }
            bufferout[col] = xx;
        }
        fwrite(&bufferout[0], sizeof(float), Ncol, SPF);
    }
    free_vector_float(bufferout);
    fclose(L1);
    fclose(PF);
    fclose(SPF);
    //fclose(txt);
    fclose(finalLocations);
    printf("\n----------Writing BMP File!----------\n");
    bmp_save(bmp,"SPF_L1(ROI).bmp");
    bmp_destroy(bmp);
    printf("\nDone!\n");
}

以及Python代码：

^{pr2}$

这是需要优化的Python代码，因为nrow和ncol的值可以并且确实达到数千个数量级。在


Tags： l1ifcolfloatnulltmpintdouble
1条回答
网友
1楼 · 发布于 2024-10-04 05:26:09

一些一般性评论：

  1. 对于python，出于多种原因，最好坚持使用PEP8。Python程序员对可读性特别挑剔，基本上都遵守社区编码准则（PEP8）。避免使用camelCase，将行保持在80列以下，去掉分号，并且可以随意地偶尔忽略这些准则，因为它们会降低可读性。

  2. 如果您使用的是numpy，这里就不需要内置的array类型。我搞不懂你为什么总是前后转换。。。

  3. 使用投影库。指定要使用的基准面和椭球体，否则坐标（东距/北距或纬度/经度）完全没有意义。

  4. 不要用一个大类作为无关事物的全部。只有几个功能没什么不好的。你不需要把它变成一个类，除非这样做是有意义的。

  5. 对numpy数组使用矢量化操作。

以下是您的性能瓶颈：

    for row in range(0,self.NRows):                                
        binvaluesL1.read(L1,self.NCols);
        binvaluesPF.read(PF,self.NCols);
        dataL1 = Num.array(binvaluesL1, dtype=Num.float);
        dataPF = Num.array(binvaluesPF, dtype=Num.float);
        dataSPF = dataL1 + dataPF;
        binvaluesSPF.fromlist(Num.array(dataSPF).tolist());                                
        for col in range(0,self.NCols):                     
            if(dataL1[col] >= self.MetalThreshold and dataPF[col] >= self.PFThreshold):
                if(col >= int(firstPoint[0]) and col <= int(secondPoint[0]) and row >= int(firstPoint[1]) and row <= int(secondPoint[1])):
                    posX = self.startX + (self.gridSizeX * col);                         
                    posY = self.startY - (self.gridSizeY * row);
                    self.UTM2LL(posY,posX);
                    tmp1 = self.decdeg2dms(posY);
                    tmp2 = self.decdeg2dms(posX);
                    strTarget = "Decimal Degree:: " + str(posX) + "E " + str(posY) + "N \t Lat long:: " + str(tmp1) + " " + str(tmp2) + "\n";
                    targetFilename.write(strTarget);
        binvaluesSPF.tofile(SPF);

你最大的问题之一是你阅读数据的方式。然后不断地把它转换成数组。完全没有必要跳过那些圈套。Numpy将像array那样为您解压二进制浮点。在

只要做grid = np.fromfile(yourfile, dtype=np.float32).reshape(ncols, nrows)。（在回路外。）

之后，只需几行代码就可以轻松地对嵌套循环进行矢量化和表示。在

下面是我如何编写你的代码。这可能不会按原样运行，因为我无法用您的数据进行测试。不过，它应该能给你一些大概的想法。在

^{pr2}$

相关问题 更多 >

    编程相关推荐