科学网

 找回密码
  注册

tag 标签: alignment

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

Rosalind - HAMM - Alignment
lucky7 2019-9-8 10:27
因为给的序列很短(要求不到1kb),我的无脑loop解法毫无压力: n=0 foriinrange(len(s1)): ifs1 !=s2 : n+=1 print(n) 当然可以用zip()函数代替for loop做迭代会更pythonic. 而且zip还有一个好处是产生iterator, 即每次内存里只放一次iteration的内容,这样在序列很长的时候就不会有内存的问题: sum( ) python3之后这些加减乘除等于不等于都可以operator模块了,确实也快: sum(map(operator.ne,s1,s2)) 再重温一下与重复相关的几个概念( 摘自 ): 循环(loop),指的是在满足条件的情况下,重复执行同一段代码。比如,while语句。 迭代(iterate),指的是按照某种顺序逐个访问列表中的每一项。比如,for语句。 遍历(traversal),指的是按照一定的规则访问树形结构中的每个节点,而且每个节点都只访问一次。 递归(recursion),指的是一个函数不断调用自身的行为。比如,以编程方式输出著名的斐波纳契数列。
个人分类: Rosalind|1673 次阅读|0 个评论
HiCPro分析流程详解
luria 2019-6-17 14:20
本篇接着上一篇《 HiCPro 的安装与使用 》,详细讲解 Hi-C 数据比对软件 HiCPro 的分析流程。 HiCPro 的安装与使用,请查看: http://blog.sciencenet.cn/blog-2970729-1182259.html 3. HiCPro 分析流程 HiCPro 处理各步骤流程如下,总体来说可以分为两大部分,对应 HiCPro 的两个 Steps : step1 比对, step2 Hi-C fragment 相关分析 3.1 HiCPro 先采用 bowtie 分别对 PE reads 进行比对 R1 reads 比对: /path/to/bowtie2 --very-sensitive -L 30 --score-min L,-0.6,-0.2 --end-to-end --reorder --un bowtie_results/bwt2_global/sample/sample_R1_sample_genome_ref.bwt2glob.unmap.fastq --rg-id BMG --rg SM:sample_R1 -p 24 -x /path/to/ref/index -U rawdata/sample/sample_R1.fastq.gz 2 logs/sample/sample_R1_bowtie2.log| /path/to/samtools view -F 4 -bS - bowtie_results/bwt2_global/sample/sample_R2_sample_genome_ref.bwt2glob.bam ## bowtie2 相关参数如下 # 注意:这里使用 --very-sensitive 。采用 bowtiew 的 Hi-C 比对软件通常会用最严格比对 # bowtie2 -L 选项是 seed substrings 的长度, 3-32 之间 # bowtie2 --end-to-end 是 entire read must align; no clipping # bowtie2 --reorder 是 force SAM output order to match order of input reads # bowtie2 --un 选项表示输出未比对上的 fastq 序列到该文件中 # bowtie2 --rg-id --rg 指定 group id 相关信息 # bowtie2 -p 指定线程数 # bowtie2 -x 指定参考基因组索引 # bowtie2 -U 指定输入的单端 reads 文件 ## samtools view 相关参数 -F 4 清理 unmapped read 。因为未比对上的 reads 已经记录在 unmap.fastq 文件中。输出的 bam 文件中不会再保留未比对上的 alignment 信息,这样可以减少后续读取文件的速度。 R2 reads 比对参数与 R1 一致 /path/to/bowtie2 --very-sensitive -L 30 --score-min L,-0.6,-0.2 --end-to-end --reorder --un bowtie_results/bwt2_global/sample/sample_R2_sample_genome_ref.bwt2glob.unmap.fastq --rg-id BMG --rg SM:sample_R2 -p 24 -x /path/to/ref/index -U rawdata/sample/sample_R2.fastq.gz 2 logs/sample/sample_R2_bowtie2.log| /path/to/samtools view -F 4 -bS - bowtie_results/bwt2_global/sample/sample_R2_sample_genome_ref.bwt2glob.bam 完成后会在 bowtie_results 目录下生成一个 bwt2_global 的文件夹,该文件夹下还有一个以样本名命名的子目录,其中包括 R1 比对 bam 文件, unmapped fastq 文件, R2 的比对 bam 文件, unmapped fastq 文件 3.2 对未比对上的 reads 进行 trim 和再比对 # 对 R1 未比对上的 reads 进行 trim /path/to/HiC-Pro_2.11.1/scripts/cutsite_trimming --fastq bowtie_results/bwt2_global/sample/sample_R1_sample_genome_ref.bwt2glob.unmap.fastq --cutsite GATCGATC --out bowtie_results/bwt2_local/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_trimmed.fastq logs/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_readsTrimming.log 21 # --fastq 指定输入的未比对上的 reads # --cutsite 指定酶切位点序列 # --out 输出 trim 之后的 reads # 对 R2 未比对上的 reads 进行 trim 。参数与 R1 一致 /path/to/HiC-Pro_2.11.1/scripts/cutsite_trimming --fastq bowtie_results/bwt2_global/sample/sample_R2_sample_genome_ref.bwt2glob.unmap.fastq --cutsite GATCGATC --out bowtie_results/bwt2_local/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_trimmed.fastq logs/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_readsTrimming.log 21 # 对 R1 trimmed reads 进行再比对 注意:这一步与第一次比对 ( 即 3.1 )节中的比对过程差异是: samtools view 没有用 -F 4 ,即将未比对上的 reads 也记录在了 bam 文件中 这是因为后期将会统计比对率,需考虑到 unmapped reads /path/to/bowtie2 --very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder --rg-id BML --rg SM:sample_R1_sample_genome_ref.bwt2glob.unmap -p 24 -x /path/to/ref/index -U bowtie_results/bwt2_local/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_trimmed.fastq 2 logs/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_bowtie2.log | /path/to/samtools view -bS - bowtie_results/bwt2_local/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_bwt2loc.bam # 对 R2 trimmed reads 进行再比对 /path/to/bowtie2 --very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder --rg-id BML --rg SM:sample_R2_sample_genome_ref.bwt2glob.unmap -p 24 -x /path/to/ref/index -U bowtie_results/bwt2_local/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_trimmed.fastq 2 logs/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_bowtie2.log | /path/to/samtools view -bS - bowtie_results/bwt2_local/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_bwt2loc.bam 运行完成后会在 bowtie_results 目录下生成一个 bwt2_local 的文件夹,该文件夹下还有一个以样本名命名的子目录,其中会生成 R1 trimmed reads 比对的 bam 结果,以及 trim 后还是无法比对上的 reads 序列文件。 R2 的也有类似的两个文件。 3.3 分别对 R1 R2 reads 两次比对的结果合并 # 对 R1 两次比对的结果合并 /path/to/samtools merge -@ 24 -n -f bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.bam bowtie_results/bwt2_global/sample/sample_R1_sample_genome_ref.0.bwt2glob.bam bowtie_results/bwt2_local/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_bwt2loc.bam # -f 选项,表示如果存在输出文件则 overwrite # -n 表示 Input files are sorted by read name # 对 R2 两次比对的结果合并 /path/to/samtools merge -@ 24 -n -f bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.bam bowtie_results/bwt2_global/sample/sample_R2_sample_genome_ref.bwt2glob.bam bowtie_results/bwt2_local/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_bwt2loc.bam # 对 R1 合并后的结果排序 /path/to/samtools sort -@ 24 -n -T tmp/sample_R1_sample_genome_ref -o bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.sorted.bam bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.bam # 对 R2 合并后的结果排序 /path/to/samtools sort -@ 24 -n -T tmp/sample_R2_sample_genome_ref -o bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.sorted.bam bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.bam # 将 R1 输出的 .sort.bam 文件改名为 .bam 文件,这一步估计是程序写完后发现有个 bug ,为了不改动后面的程序,加的一步 mv bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.sorted.bam bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.bam # 将 R2 输出的 .sort.bam 文件改名为 .bam 文件 mv bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.sorted.bam bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.bam 运行完成后,会在 bowtie_results 目录下生成一个 bwt2 文件夹,和以样本名为名字的子文件夹。在其中生成两个以 bwt2merged.bam 结尾的文件,分别代表 R1 R2 的结果 3.4 统计比对率 采用的方法是 samtools view -c ,具体如下: # 分别统计 R1 R2 最终结果中总的 reads 和成对的 reads /path/to/samtools view -c bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.bam /path/to/samtools view -c bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.bam /path/to/samtools view -c -F 4 bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.bam /path/to/samtools view -c -F 4 bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.bam # 分别统计 R1 R2 两次比对结果中,比对上的 reads 数 /path/to/samtools view -c -F 4 bowtie_results/bwt2_global/sample/sample_R1_sample_genome_ref.bwt2glob.bam /path/to/samtools view -c -F 4 bowtie_results/bwt2_global/sample/sample_R2_sample_genome_ref.bwt2glob.bam /path/to/samtools view -c -F 4 bowtie_results/bwt2_local/sample/sample_R1_sample_genome_ref.bwt2glob.unmap_bwt2loc.bam /path/to/samtools view -c -F 4 bowtie_results/bwt2_local/sample/sample_R2_sample_genome_ref.bwt2glob.unmap_bwt2loc.bam 3.5 采用 HiCPro 的 mergeSAM.py 程序合并 PE reads /path/to/python /path/to/HiC-Pro_2.11.1/scripts/mergeSAM.py -q 10 -t -v -f bowtie_results/bwt2/sample/sample_R1_sample_genome_ref.bwt2merged.bam -r bowtie_results/bwt2/sample/sample_R2_sample_genome_ref.bwt2merged.bam -o bowtie_results/bwt2/sample/sample_sample_genome_ref.bwt2pairs.bam # -q 指定最小的 mapping quality # -t 表示需生成比对统计结果 # -v 用于 debug # -f/--forward -r/--reverse 分别输入上一步生成的 PE reads 分开比对最终的 bam 文件 ( 即 3.3 中结果 ) # -o 输出处理后的结果 此外, HiCPro 的 mergeSAM.py 程序还有以下几个参数: report singleton 给出单端比对的结果 report multiple hits 给出多处比对的结果 运行完成后,会在 bwt2 的子目录中生成后缀为 bwt2pairs.bam 的过滤结果文件。以及其统计文件(后缀为 .bwt2pairs.pairstat 的文件),其中会统计完整的比对信息。包括 Total_pairs_processed, Unmapped_pairs, Low_qual_pairs, Unique_paired_alignments, Multiple_pairs_alignments, Pairs_with_singleton, Low_qual_singleton, Unique_singleton_alignments, Multiple_singleton_alignments, Reported_pairs 至此,比对的部分全部完成,所有的结果都在 bowtie_results 目录下 ===================================================== 再进行 Hi-C fragment 相关分析,所有的结果都在 hic_results 目录下 3.6 利用 HiCPro 的 mapped_2hic_fragments.py 程序将比对结果转化为 Hi-C 片段信息 /path/to/python /path/to/HiC-Pro_2.11.1/scripts/mapped_2hic_fragments.py -v -a -f /path/to/restriction_enzyme_cutting_site.MboI.txt -r bowtie_results/bwt2/sample/sample_sample_genome_ref.bwt2pairs.bam -o hic_results/data/sample # -v 用于 debug # -a 记录所有的信息,包括 self-circle, dangling end 等 # -f 指定最开始时,检测出的基因组序列上的酶切位点信息文件,即 HiC-Pro_2.11.1/bin/utils/digest_genome.py 生成的结果 # -r 指定 bowtie2 比对的最终结果,即 3.5 节中的结果 # -o 指定输出目录 此处,还可以指定 insert size 的阈值等,具体可参见 --help 再对输出的 valid pairs 文件进行排序: LANG=en; sort -T tmp -k2,2V -k3,3n -k5,5V -k6,6n -o hic_results/data/sample/sample_sample_genome_ref.bwt2pairs.validPairs hic_results/data/sample/sample_sample_genome_ref.bwt2pairs.validPairs # 这里是直接将 valid paris 原路径排序 3.7 对所有的 valid pairs 进行合并,并且去掉 PCR duplication LANG=en; sort -T tmp -S 50% -k2,2V -k3,3n -k5,5V -k6,6n -m hic_results/data//sample/sample_sample_genome_ref.bwt2pairs.validPairs | awk -F\\t 'BEGIN{c1=0;c2=0;s1=0;s2=0}(c1!=$2 || c2!=$5 || s1!=$3 || s2!=$6){print;c1=$2;c2=$5;s1=$3;s2=$6}' hic_results/data/sample/sample.allValidPairs 这一步是先将多个 valid pairs 文件进行合并,例如加测了几次,如果之前没有合并,这里可以合到一起。然后再确定当前行和上一行是否相同,如果相同则为 PCR duplication, 需去除 3.6 和 3.7 节的结果都生成在 hic_results/data 目录中 3.8 采用 HiCPro 的 merge_statfiles.py 程序对 bowtie2 比对的多个统计结果合并 /path/to/python /path/to/HiC-Pro_2.11.1/scripts/merge_statfiles.py -d bowtie_results/bwt2/sample/ -p *_R1*.mapstat -v hic_results/stats/sample/sample_R1.mmapstat /path/to/python /path/to/HiC-Pro_2.11.1/scripts/merge_statfiles.py -d bowtie_results/bwt2/sample/ -p *_R2*.mapstat -v hic_results/stats/sample/sample_R2.mmapstat /path/to/python /path/to/HiC-Pro_2.11.1/scripts/merge_statfiles.py -d bowtie_results/bwt2/sample/ -p *.pairstat -v hic_results/stats/sample/sample.mpairstat /path/to/python /path/to/HiC-Pro_2.11.1/scripts/merge_statfiles.py -d hic_results/data//sample/ -p *.RSstat -v hic_results/stats/sample/sample.mRSstat # 这一步统计的结果都在 hic_results/stats 目录下 3.9 跟据 BIN_SIZE 来构建 matrix 这一步会按照 HiCPro config-hicpro.txt 文件中指定的 BIN_SIZE ,对 valid pairs 进行分配,构建 matrix cat hic_results/data//sample/sample.allValidPairs | /path/to/HiC-Pro_2.11.1/scripts/build_matrix --matrix-format upper --binsize 20000 --chrsizes /path/to/ref/reference.size --ifile /dev/stdin --oprefix hic_results/matrix/sample/raw/20000/sample_${bsize} # --chrsizes 是最初统计出的基因组中每条序列的长度 运行完,各分辨率的 matrix 输出到 hic_results/matrix/sample/raw 目录下 3.10 对统计结果画图 会生成 5 个图,分别是 plotHiCContactRanges_sample.pdf, plotHiCFragmentSize_sample.pdf, plotMappingPairing_sample.pdfplot,HiCFragment_sample.pdf, plotMapping_chrysanthemum.pdf, 如下图: 3.11 采用 ice 对 raw matrix 做 normalization /path/to/python /path/to/HiC-Pro_2.11.1/scripts/ice --results_filename hic_results/matrix/sample/iced/20000/sample_20000_iced.matrix --filter_low_counts_perc 0.02 --filter_high_counts_perc 0 --max_iter 100--eps 0.1 --remove-all-zeros-loci --output-bias 1 --verbose 1 hic_results/matrix//sample/raw/20000/sample_20000.matrix # 对各分辨率生成的 raw matrix 做 normalization ,结果都输出在 hic_results/matrix/sample/ice 目录下
个人分类: Hi-C|13700 次阅读|0 个评论
[转载]Kart -- 一个分而治之NGS read比对算法,适合于long reads
chinapubmed 2019-3-14 08:29
下一代测序(Next-generation sequencing,NGS)提供了巨大的机遇以核苷酸分辨率研究全基因组变异。由于巨大的数据量,NGS应用需要非常快速且准确的比对算法。大部分现有read比对算法基本上都采用种子和扩展策略,这本质上是顺序的,并且在更长reads上耗费更多时间。 我们开发了一种分而治之(divide-and-conquer)算法,叫做Kart,它通过将一条read分割成能够独立比对的小片段能够像短reads一样快地处理长reads。我们的实验结果表明更耗时的空位比对需要的平均片段大小为约20bp,而不管原始read长度。此外,它能够容忍更高的错误率。实验显示当错误率高达15%时,Kart比其他比对器在更长reads上花费更少的时间,并且仍然产生可靠的比对。 ABCD为完全匹配的simple pair 下载:https://github.com/hsinnan75/Kart/ 语言:C++ 时间:20170405 参考: Kart: a divide-and-conquer algorithm for NGS read alignment
个人分类: 软件|2077 次阅读|0 个评论
序列多重比对工具:MUSCLE
vesperlight 2017-8-15 09:27
Muscle MUSCLE是RC Edgar开发的序列多重比对(Multiple Sequence Alignment,MSA)工具 下载和相关说明地址为http://www.drive5.com/muscle/manual/ 1、比对并保存比对结果为Fasta格式文件 muscle -in seqs.fa -out seqs.afa 对于大数据集可以使用 muscle -in seqs.fa -out seqs.afa -maxiters 2 MUSCLE默认使用高准确度的比对方式,若需要更快但精度较低的方法可以使用: 用于氨基酸序列比对 muscle -in seqs.fa -out seqs.afa -maxiters 1 -diags -sv -distance1 kbit20_3 用于核算序列比对 muscle -in seqs.fa -out seqs.afa -maxiters 1 -diags 2、将比对结果转换为CLUSTALW格式的文件 muscle -in seqs.fa -clw 参数 -in 输入文件必须为fasta格式,如果序列中存在gaps,gaps将会被丢弃 -out 输出文件 3、根据多重比对结果构建UPGMA树 muscle -maketree -in seqs.afa -out seqs.phy 4、根据多重比对结果构建NJ(Neighbor-Joining)树 muscle -maketree -in seqs.afa -out seqs.phy -cluster neighborjoining 5、对已有比对结果加入新的序列 已有一个msa结果,想加入一条新的序列 muscle -profile -in1 existing_msa.afa -in2 new_seq.fa -out combined.afa 6、如果序列有多条,则先对需要加入的序列进行多重比对,然后对两个多重比对结果进行比对(同下) muscle -in new_seqs.fa -out new_seqs.afa muscle -profile -in1 existing_aln.afa -in2 new_seqs.afa -out combined.afas 7、两个比对结果进行比对 muscle -profile -in1 one.afa -in2 two.afa -out both.afa 8、提炼已有MSA muscle -in msa.afa -out refined_msa.afa -refine
个人分类: 生信|20165 次阅读|0 个评论
基于HOG特征的SDM face alignment
wanglin193 2017-4-17 18:06
测试了基于梯度下降方法,SDM是根据当前形状位置提取的图像信息,直接用回归法计算出形状更新增量的方法,即dShape_i = M_i*dFeature_i,其中M_i是计算第i次回归时使用的矩阵。 很大程度上,原始SDM算法中HOG特征起了很大作用。用原始图像patch的像素值直接来做回归效果并不好,而基于梯度方向直方图的特征具有光照无关的特性,类似Garbor小波,是具有很高冗余度的信息,即特征矢量的长度往往大于patch的像素值。 68个脸部特征点的HOG串联的矢量作为输入,左乘M矩阵,输出136维(68×2)的xy图像位置增量,M矩阵通过给定已知输入和输出的样本集得到。求解M的最小二乘解可以比较耗时,所以也可以先把样本集特征数据做个PCA降维,至少SDM论文里是这样写的。 源码放在这里: https://github.com/wanglin193/SupervisedDescentMethod HOG/SIFT等手工选取的特征应进一步去除冗余,可以用训练的方法选取最适合的特征矢量,比如LBF,或者深度学习方法提取的特征。
个人分类: ASM/AAM|4807 次阅读|0 个评论
PSL format for blat software
ljxue 2014-11-10 02:56
https://www.biostars.org/p/44633/ In general the coordinates in psl files are “zero based half open.” The first base in a sequence is numbered zero rather than one. When representing a range the end coordinate is not included in the range. Thus the first 100 bases of a sequence are represented as 0-100, and the second 100 bases are represented as 100-200. There is a another little unusual feature in the .psl format. It has to do with how coordinates are handled on the negative strand. In the qStart/qEnd fields the coordinates are where it matches from the point of view of the forward strand (even when the match is on the reverse strand). However on the qStarts[] list, the coordinates are reversed. Here's an example of a 30-mer that has 2 blocks that align on the minus strand and 2 blocks on the plus strand (this sort of stuff happens in real life in response to assembly errors sometimes). 0 1 2 3 tens position in query 0123456789012345678901234567890 ones position in query ++++ +++++ plus strand alignment on query -------- ---------- minus strand alignment on query Plus strand: qStart 12 qEnd 31 blockSizes 4,5 qStarts 12,26 Minus strand: qStart 4 qEnd 26 blockSizes 10,8 qStarts 5,19 Essentially the minus strand blockSizes and qStarts are what you would get if you reverse complemented the query. However the qStart and qEnd are non-reversed. To get from one to the other: qStart = qSize - revQEnd qEnd = qSize - revQStart
个人分类: Bioinformatics|3743 次阅读|0 个评论
ASM/AAM手工标定的Matlab脚本
热度 9 wanglin193 2011-5-25 18:57
ASM/AAM手工标定的Matlab脚本
有段时间我一直满足于在脸上画个框,后来又满足于在脸上标些点,再后来在这方面就没有进展了。剩下不多的东西包括硬盘某个文件夹里不同时期下载的paper:光流、特征检测、ASM、AAM...... 还有各种版本的人脸检测和跟踪的代码,里面藏着很多垃圾,比如每个版本的文件夹都包含一个debug和release目录。相对来说那些Matlab小程序就挺不错,不太占地方,用几行代码就能实现挺复杂的功能,比如下面这个AAM或者ASM训练用的手工标点的脚本,我把它放在这里晒晒。 程序下载: Edit68CMU_pack.rar 对于人脸检测的人脸样本准备,通常需要手工从图像中截取人脸区域,为了程序鲁棒性的考虑,生成的人脸样本需要还要有一定角度的旋转;对基于统计训练的图像自动标定(image alignment,比如Active shape model 或Active appearance model)方法来说,第一步工作就是手工在图上把物体轮廓标出来,这些点一般称为landmark。在没有用过OpenCV之前,要写一个支持鼠标输入交互选点的程序还是挺麻烦的,但是等到我开始使用Matlab,并且知道有个ginput函数后,这事就变得无比简单了。 首先,因为ASM的形状(shape)点是有序号的,每个点都对应一个目标上的关键点,是有具体含义的。一个标定的方法是用sh=ginput(68)一口气选取68个点,但是这中间一旦出错就前功尽弃了。所以比较好的办法是先把一个差不多的标准形状铺在输入的图像上,然后交互地调整每个点的位置,直到每个点都放到对应的位置为止,这样就把"标定形状"变成"编辑形状"了。其次,如果一个人脸图像"A.jpg"对应的landmark已经标记好了,并且保存在一个比如叫做"A.lmk"的文件中,过了一阵子我又想重新修改某些点的位置,或者因为不同的应用,形状的配置发生变化需要增加删除某些点,程序也应该满足编辑lmk文件个功能。 下面的程序就主要完成这两个功能: 1.先准备个标准形状做参考,比如从论文里找这样一张图,用ginput按顺序选68个点(68个点的定义见CMU那帮人的论文,文献 ,也可以用画图软件事先标好你要的点,比如下图的蓝点,这样的好处是可以保证使用ginput(68)一次成功): imshow('meanface.jpg'); sh = ginput(68); 这样这些点的X,Y坐标就放在sh的两列里。把它作为一个标准形状,存在一个文本文件里。 2.程序开始,先要load标准形状。然后通过对话框选取一个图像文件,并显示 = uigetfile( , 'Pick an image'); Aim = imread( ); imshow(Aim, 投影到图像中。 4.如果对应的lmk文件已经存在,则读入lmk,把形状 铺在图像上。 5.根据上面的输入信息,粗略计算一个人脸区域box,并且放大显示 im=imcrop(Aim,box); imshow(im, =ginput(1); %最近点的序号存在ind里 =min((xc-Xo).^2 + (yc-Yo).^2); 5.2 再用ginput(1)选择要移动的位置,替换旧的值 = ginput(1); Xo(ind)=xc; Yo(ind)=yc; 重复这两个步骤,直到满意为止。 6. 保存形状文件(.lmk) fid = fopen( , 'w'); fprintf(fid,'%6.2f %6.2f\r\n', '); fclose(fid); 另外利用delaunay函数事先计算平均形状的三角网格(mesh): mtri=delaunay(meanshape68(:,1),meanshape68(:,2)); mtri的每行存储的是3角形的3个顶点序号,我们也可以利用这个三角剖分信息,把网格画在人脸上。 遗憾的是这3个点不都是顺时针或逆时针的,也就是说delaunay输出的这些三角形面片的法线方向不是朝一个方向的,这就需要交换某些三角形内某对点的位置,网格三角剖分工作甚至可以不用delaunay函数,而用手工完成,因为这个步骤都是一次性的。这些三角网格的主要用处是把人脸图像纹理warp到平均形状,叫做shape-free patch,这是AAM的关键步骤之一。如下图所示。左边是原始图的网格,右上角是表示在平均形状位置的每个三角形有一组仿射变换参数(用来到原始图像中采样,每个三角形内的图像采样点相对于该三角形网格的参数是固定的,也只要计算一次即可),右下角是得到的shape-free patch。 这样做,可以有助于实现形状(shape)和纹理(appearance)的的分别训练,具体更进一步实现自动人脸校准(face alignment),就需要好好看看下面两个文章了。 : 1. AAM/ASM创始人Tim Cootes的文章: http://personalpages.manchester.ac.uk/staff/timothy.f.cootes/refs.html 2.CMU的方法,以Lucas-Kanade 算法为基础(同时可参考我之前贡献的Lucas-Kanade 算法的文档) Active Appearance Models Revisited : http://www.ri.cmu.edu/publication_view.html?pub_id=4601 3.3D电影Avatar就是用了类似的方法进行人脸表情跟踪的。 http://www.meraforum.com/showthread.php?t=36119
个人分类: ASM/AAM|26191 次阅读|16 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-5-19 19:36

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部