脉冲式核磁数据中计算基线高度：在大量点中找到基线高度，C#实现

ゝ一纸荒年。 2022-11-29 05:43 18阅读 0赞

做化学好多年，一直想为Researchers做点贡献。一直梦想着写个核磁软件，无奈无从下手。最近看到http://blog.sina.com.cn/s/blog\_1349667240102wzrk.html上说起核磁文件的格式，我才有了开始写核磁软件的可能。写了两天了，能实现显示、积分的功能了。现在要寻找谱图的基线。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dhbmdlaWwwMDc_size_16_color_FFFFFF_t_70_pic_center]  
核磁数据中1r中储存了65536个Y轴数据。其中X轴的数据是在-4~16之间均匀分布的。把所有这些数据画在一张图上就得到了上述核磁的原始谱图。怎样才能得到这个基线的高度呢。基线高度肯定不是谱图中最低的点的高度，也可能不是最高点的高度。 (屏幕左上角是坐标原点，Y值越大，越是在下面)。

可以确定的一点是基线，或者基线附近的点是最多的，真正产生峰的点是比较少的。所以，我们可以根据高度进行二分：我们先把所有原始数据放到listPoint中。  
1、找到listPoint数据中的Y的最大值和最小值。  
2、根据Y的最大值和最小值，把Y值分成高度相等的2个区域，  
3、遍历listPoint数据中的Y值，把数据放到指定的区域中。清空listPoint中的数据。  
4、找到数据比较多的区域，并把这个区域中的点放到listPoint中。  
5、重复1，2，3，4，直到数据比较多的区域中的最大值和最小值相等。基线高度就是这个区域中的最大值 (最大值、最小值、平均值都相等)。  
每次遍历的时间复杂度为O(N)。如果谱图中任意两个点之间的高度的最大差为M，用二分的方法缩减高度差的值，一共需要O(logM)次，所以整个过程的时间复杂度为O(N \* logM)。这个寻找基线的方法的效率比较高，其实我也想不到其它方法找基线了 。这个方法可能是有缺陷的：如果基线很不平整，这个方法就可能有问题。例如样品中有一些亚铁离子等顺磁性的物质，核磁结果就很难看，其中一个表现就是基线很不平整。。。，另外，我们的假设是基线附近的点最多。但是如果基线附近的点不是最多，是不是也会出现异常？  
值得注意的一点：二分过程中，随着Y变化，每次需要遍历的元素个数也是在不断减小的，但是遍历时间复杂度不变。从下图可以看到，每进行一次二分，总归有一些点被丢弃，而且每次丢弃的点的个数是不断增加的。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dhbmdlaWwwMDc_size_16_color_FFFFFF_t_70_pic_center 1]  
下面是实现代码：

public static int getBaseLineY()
    {
        
        int pointNum = R.listPoint.Count; //点的个数
        List<Point> listPoint = new List<Point>();
        listPoint.AddRange(R.listPoint); //避免用R.listPoint
    
        int minY = 9999999;
        int maxY = 0;
    
        //找到minY、maxY
        for (int i = 0; i < pointNum; i++)
        {
        
            if (maxY < listPoint[i].Y)
            {
        
                maxY = listPoint[i].Y;
                //continue;
            }
    
            if (minY > listPoint[i].Y)
            {
        
                minY = listPoint[i].Y;
            }
        }
    
        List<Point> listPointMin = new List<Point>(); //存放比中间Y值小的点
        List<Point> listPointMax = new List<Point>(); //存放比中间Y值大的点
    
        //开始二分计算
        int minY1 = minY; //第一部分的最大值和最小值
        int maxY1 = maxY;
    
        int minY2 = minY; //第二部分的最大值和最小值
        int maxY2 = maxY;
    
        int middle = 0; //中间值
        //while ((maxY1 - minY1) * 200 > (maxY - minY)) //控制精度的方式不好。数据有可能越界
        while (maxY1 > minY1)//一共需要循环logN次
        {
        
            middle = minY1 + (maxY1 - minY1) / 2;//找到中间值
            pointNum = listPoint.Count; //更新pointNum
    
            //Console.WriteLine("maxY1 = " + maxY1 + ", minY1 = " + minY1 + ", maxY1 - minY1 = " + (maxY1 - minY1) + ", maxY - minY = " + (maxY - minY) + ", middle = " + middle + ", pointNum = " + pointNum + ", listPointMax = " + listPointMax.Count + ", listPointMin = " + listPointMin.Count);
            listPointMax.Clear();
            listPointMin.Clear();
    
            for (int i = 0; i < pointNum; i++) //循环N次
            {
        
                if (middle < listPoint[i].Y) //按照Y的值进行分组
                {
        
                    listPointMax.Add(listPoint[i]);
    
                    if (maxY2 < listPoint[i].Y)
                    {
        
                        maxY2 = listPoint[i].Y;
                    }
    
                    else if (minY2 > listPoint[i].Y)
                    {
        
                        minY2 = listPoint[i].Y;
                    }
    
                }
                else //(middle >= listPoint[i].Y)
                {
        
                    listPointMin.Add(listPoint[i]);
    
                    if (maxY1 < listPoint[i].Y)
                    {
        
                        maxY1 = listPoint[i].Y;
                    }
    
                    else if (minY1 > listPoint[i].Y)
                    {
        
                        minY1 = listPoint[i].Y;
                    }
                }
            }
    
            listPoint.Clear();
    
            //完成一次统计，重置listPoint为数量比较多的那个listPoint
            if (listPointMax.Count > listPointMin.Count)
            {
        
                listPoint.AddRange(listPointMax);
                minY1 = middle + 1; //二分的关键
                maxY1 = maxY2;
            }
            else
            {
        
                listPoint.AddRange(listPointMin);
                minY2 = minY1; //二分的关键
                maxY2 = middle - 1;
            }
        }
        return maxY1;
    }

运行结果如下：

![在这里插入图片描述][20200819143406816.png_pic_center]正如上面所说的，每次都会丢弃一些点，且每次丢弃的点的数量越来越多。49->101->408->869->905->706(这个数据小于905，所以证明丢弃的点越来越多这句话是错的)->2093。

高度差从338开始减半衰减，只用了8次就到0，得到基线高度为600。最后，利用这个返回值，在图上画出一条绿色的直线 (如最上面的图片所示)。

这个圈子绕的够大的！！！

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dhbmdlaWwwMDc_size_16_color_FFFFFF_t_70_pic_center]: /images/20221124/b45b07fe88964060a3554039a027848b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dhbmdlaWwwMDc_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221124/93eb98640d9e42e8b8b42939ddc4147a.png
[20200819143406816.png_pic_center]: /images/20221124/0e8ecd23192b40f9bafa93c493a907e9.png