LAB 07 · 空间格、方向桶与传统表示

import json
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=4, suppress=True)
plt.rcParams.update({'figure.dpi':110,'font.size':10,'axes.spines.top':False,'axes.spines.right':False})
BLUE, GOLD, GREEN, PURPLE = '#002676','#FDB515','#00856A','#69458D'
CHECKS={}; RESULTS={}
def check(name,condition):
    CHECKS[name]=bool(condition)
    assert CHECKS[name], name
def table(headers,rows,title):
    fig,ax=plt.subplots(figsize=(max(7,len(headers)*1.5),1+len(rows)*.36))
    ax.axis('off'); ax.set_title(title,pad=12,color=BLUE)
    tab=ax.table(cellText=rows,colLabels=headers,loc='center',cellLoc='center')
    tab.auto_set_font_size(False);tab.set_fontsize(10);tab.scale(1,1.45)
    for (r,c),cell in tab.get_celld().items():
        cell.set_edgecolor('#CBD5DF')
        if r==0:cell.set_facecolor('#E9F0F7')
    fig.tight_layout();plt.show()

P=dict(sift_rotation=0.0, dominant_orientation=0.0,
       sift_magnitudes=[1.0,2.0,1.0,2.0], neighbor_cell=[2,4],
       hog_angles=[10.0,190.0,90.0],hog_magnitudes=[2.0,2.0,1.0],
       hog_offset=0.0,hog_gain=1.0,normalization_epsilon=1e-6,
       pyramid_grid=[2,2])
P.update(globals().get('USER_OVERRIDES',{}))
print('Teaching configuration:',json.dumps(P,ensure_ascii=False))

centers=np.arange(8)*45.
angles=(np.array([0.,45.,45.,90.])+float(P['sift_rotation'])-float(P['dominant_orientation']))%360
mags=np.asarray(P['sift_magnitudes'],dtype=float);assert mags.shape==(4,) and np.all(mags>=0)
def hard_signed_vote(theta,magnitude):
    ids=np.floor((np.asarray(theta)%360)/45.+.5).astype(int)%8
    return np.bincount(ids,weights=magnitude,minlength=8)
sift_hist=hard_signed_vote(angles,mags)
check('SIFT_vote_mass',np.isclose(sift_hist.sum(),mags.sum()))
if np.array_equal(angles,[0,45,45,90]) and np.array_equal(mags,[1,2,1,2]):check('SIFT_slide_r2c3',np.array_equal(sift_hist,[1,3,2,0,0,0,0,0]))
fig,axes=plt.subplots(1,2,figsize=(10,4));origins=np.array([[.3,.35],[1.45,.35],[.3,1.45],[1.45,1.45]])
for (x,y),a,m in zip(origins,angles,mags):
    rad=np.deg2rad(a);axes[0].annotate('',xy=(x+.65*np.cos(rad),y+.65*np.sin(rad)),xytext=(x,y),arrowprops=dict(arrowstyle='-|>',color=BLUE,lw=2))
    axes[0].text(x-.1,y-.25,f'{a:g} deg; m={m:g}',fontsize=9)
axes[0].set(xlim=(-.2,2.4),ylim=(-.2,2.4),aspect='equal',title='Constructed gradients in cell r2c3',xlabel='x (right)',ylabel='y (up)');axes[0].grid(alpha=.2)
axes[1].bar(centers,sift_hist,width=30,color=BLUE);axes[1].set(xticks=centers,xlabel='Relative direction (degrees)',ylabel='Sum of magnitudes',title='Eight signed direction bins')
fig.tight_layout();plt.show()
table(['bin b (zero-based)','direction','vote'],[[i,f'{a:g} deg',f'{v:g}'] for i,(a,v) in enumerate(zip(centers,sift_hist))],'Local direction histogram')
RESULTS['sift_vote']={'relative_angles':angles.tolist(),'histogram':sift_hist.tolist()}

def address(r,c,b):
    assert 1<=r<=4 and 1<=c<=4 and 0<=b<8
    return ((r-1)*4+c-1)*8+b
nr,nc=map(int,P['neighbor_cell']);assert 1<=nr<=4 and 1<=nc<=4 and (nr,nc)!=(2,3)
neighbor_hist=hard_signed_vote([45,225],[1,2])
descriptor=np.zeros((4,4,8));descriptor[1,2]=sift_hist;descriptor[nr-1,nc-1]=neighbor_hist
flat=descriptor.ravel();focus_j=address(2,3,1);neighbor_j=address(nr,nc,1)
check('SIFT_49_address',focus_j==49)
check('SIFT_neighbor_address_distinct',neighbor_j!=focus_j)
check('SIFT_flatten_roundtrip',np.array_equal(flat.reshape(4,4,8),descriptor))
check('SIFT_128_dimensions',len(flat)==128)
fig,axes=plt.subplots(1,2,figsize=(11,4))
for r in range(4):
    for c in range(4):
        fill='#FFF0B7' if (r,c)==(1,2) else '#DDECE7' if (r,c)==(nr-1,nc-1) else '#F3F6FA'
        axes[0].add_patch(plt.Rectangle((c,3-r),1,1,facecolor=fill,edgecolor='#AABBCD'))
        axes[0].text(c+.5,3-r+.66,f'r{r+1}c{c+1}',ha='center');axes[0].text(c+.5,3-r+.28,f'{address(r+1,c+1,0)}-{address(r+1,c+1,7)}',ha='center',fontsize=9)
axes[0].set(xlim=(0,4),ylim=(0,4),aspect='equal',title='4 x 4 cells; each stores eight bins');axes[0].axis('off')
axes[1].bar(np.arange(128),flat,color=BLUE,width=.8);axes[1].axvline(focus_j,color=GOLD,lw=2,label='r2c3, b1 -> 49');axes[1].axvline(neighbor_j,color=GREEN,lw=2,label=f'neighbor b1 -> {neighbor_j}')
axes[1].set(xlim=(0,127),xlabel='Descriptor index (zero-based)',ylabel='Vote',title='Same direction, different spatial addresses');axes[1].legend(fontsize=8);fig.tight_layout();plt.show()
table(['cell','bin','index','value'],[['r2c3',1,focus_j,float(flat[focus_j])],[f'r{nr}c{nc}',1,neighbor_j,float(flat[neighbor_j])]],'Address and value are different concepts')
RESULTS['sift_address']={'index49_value':float(flat[49]),'neighbor_index':neighbor_j,'descriptor_size':int(flat.size)}

def unsigned_vote(theta,magnitude,bins=9):
    theta=np.asarray(theta,dtype=float)%180.;magnitude=np.asarray(magnitude,dtype=float)
    assert theta.shape==magnitude.shape and np.all(magnitude>=0)
    q=theta/(180./bins);lo=np.floor(q).astype(int)%bins;fraction=q-np.floor(q)
    hist=np.zeros(bins);np.add.at(hist,lo,magnitude*(1-fraction));np.add.at(hist,(lo+1)%bins,magnitude*fraction)
    return hist
ha=np.asarray(P['hog_angles'],dtype=float)+float(P['hog_offset']);hm=np.asarray(P['hog_magnitudes'],dtype=float)
hog=unsigned_vote(ha,hm);rotated=unsigned_vote(ha+180,hm)
gain=float(P['hog_gain']);assert gain>0;eps=float(P['normalization_epsilon']);assert eps>0
scaled=hog*gain
normalize=lambda a:a/np.sqrt(np.dot(a,a)+eps**2)
check('HOG_unsigned_180_invariance',np.allclose(hog,rotated))
check('HOG_vote_mass',np.isclose(hog.sum(),hm.sum()))
check('HOG_10deg_two_votes',np.allclose(unsigned_vote([10],[2]),[1,1,0,0,0,0,0,0,0]))
fig,axes=plt.subplots(1,2,figsize=(11,4));bins=np.arange(9)*20
axes[0].bar(bins-3,hog,6,color=BLUE,label='base');axes[0].bar(bins+3,scaled,6,color=GOLD,label='gain-scaled')
axes[1].bar(bins-3,normalize(hog),6,color=BLUE,label='base normalized');axes[1].bar(bins+3,normalize(scaled),6,color=GREEN,label='scaled normalized')
for ax in axes:ax.set(xticks=bins,xlabel='Unsigned bin center (degrees)',ylabel='Vote');ax.legend(fontsize=8)
axes[0].set_title('Before normalization');axes[1].set_title('One-vector L2 mechanism');fig.tight_layout();plt.show()
table(['input angle','angle mod 180','magnitude'],[[f'{a:g}',f'{a%180:g}',f'{m:g}'] for a,m in zip(ha,hm)],'Direction polarity is discarded')
print('Raw votes:',hog,' normalized:',normalize(hog))
print('Classic window count: 64 x 128 px / 8 px cells -> 8 x 16 cells; 2 x 2-cell blocks -> 7 x 15 x 36 =',7*15*36)
RESULTS['hog']={'histogram':hog.tolist(),'normalized':normalize(hog).tolist(),'unsigned_180_error':float(np.max(np.abs(hog-rotated))),'illustrative_window_dimensions':3780}

training_groups=[np.array([[0,0],[0,2],[2,0],[2,2]]),np.array([[4,3],[4,5],[6,3],[6,5]])]
codebook=np.stack([g.mean(0) for g in training_groups])
new_descriptors=np.array([[2.,1.],[4.,4.]]) # 学生可改这两个新局部描述；不要重建词典。
distance=((new_descriptors[:,None,:]-codebook[None,:,:])**2).sum(2)
assignment=distance.argmin(1);counts=np.bincount(assignment,minlength=len(codebook))
table(['new descriptor','nearest center ID','squared distances'],[[str(d.tolist()),int(a)+1,str(ds.tolist())] for d,a,ds in zip(new_descriptors,assignment,distance)],'Training centers retained: (1,1), (5,4)')
SLIDE_WORDS_A=np.array(['A','A','A','B','A','A','B','C','B','B','C','C','B','C','A','B']).reshape(4,4)
WORDS_A=SLIDE_WORDS_A.copy() # 可在这里改位置或组成；例如 WORDS_A[0,0]='C'。
assert WORDS_A.shape==(4,4) and np.isin(WORDS_A,['A','B','C']).all()
WORDS_B=np.vstack([WORDS_A[2:],WORDS_A[:2]])
vocabulary=['A','B','C']
count_words=lambda a:np.array([(a==k).sum() for k in vocabulary])
bag_a,bag_b=count_words(WORDS_A),count_words(WORDS_B)
check('row_permutation_preserves_bag',np.array_equal(bag_a,bag_b))
print('Two layouts actually differ:',not np.array_equal(WORDS_A,WORDS_B))
if np.array_equal(WORDS_A,SLIDE_WORDS_A):check('slide_bag_counts',np.array_equal(bag_a,[6,6,4]))
fig,axes=plt.subplots(1,3,figsize=(11,3.5));colors={'A':BLUE,'B':GOLD,'C':GREEN}
for ax,words,title in zip(axes[:2],[WORDS_A,WORDS_B],['Layout A','Layout B: swap row halves']):
    for r in range(4):
        for c in range(4):
            word=words[r,c];ax.add_patch(plt.Rectangle((c,3-r),1,1,facecolor=colors[word],edgecolor='white'));ax.text(c+.5,3-r+.5,word,ha='center',va='center',color='black' if word=='B' else 'white',fontsize=14)
    ax.set(xlim=(0,4),ylim=(0,4),aspect='equal',title=title);ax.axis('off')
axes[2].bar(np.arange(3)-.17,bag_a,.34,color=BLUE,label='layout A');axes[2].bar(np.arange(3)+.17,bag_b,.34,color=GOLD,label='layout B');axes[2].set(xticks=np.arange(3),xticklabels=vocabulary,ylabel='Count',title='Same global bag');axes[2].legend();fig.tight_layout();plt.show()
RESULTS['bag']={'centers':codebook.tolist(),'assignments':assignment.tolist(),'new_counts':counts.tolist(),'layout_A_counts':bag_a.tolist(),'layout_B_counts':bag_b.tolist()}
# 真正的 Lloyd fit 已独立迁入第三讲 CLS02 前史选读；本节只保留预设均值与布局计数。

upper_a,lower_a=count_words(WORDS_A[:2]),count_words(WORDS_A[2:])
upper_b,lower_b=count_words(WORDS_B[:2]),count_words(WORDS_B[2:])
if np.array_equal(WORDS_A,SLIDE_WORDS_A):check('slide_upper_lower_counts',np.array_equal(upper_a,[5,2,1]) and np.array_equal(lower_a,[1,4,3]))
check('swapped_regions',np.array_equal(upper_a,lower_b) and np.array_equal(lower_a,upper_b))
table(['layout','global A/B/C','upper A/B/C','lower A/B/C'],[['A',str(bag_a.tolist()),str(upper_a.tolist()),str(lower_a.tolist())],['B',str(bag_b.tolist()),str(upper_b.tolist()),str(lower_b.tolist())]],'Coarse spatial addresses restore a distinction')
gr,gc=map(int,P['pyramid_grid']);assert gr>0 and gc>0 and 4%gr==0 and 4%gc==0
def pyramid(words,rows,cols):
    result=[count_words(words)];names=['global']
    for r in range(rows):
        for c in range(cols):
            patch=words[r*4//rows:(r+1)*4//rows,c*4//cols:(c+1)*4//cols]
            result.append(count_words(patch));names.append(f'r{r+1}c{c+1}')
    return np.concatenate(result),names
pa,regions=pyramid(WORDS_A,gr,gc);pb,_=pyramid(WORDS_B,gr,gc)
check('pyramid_dimensions',len(pa)==3*(1+gr*gc))
check('pyramid_count_mass',pa.sum()==pb.sum()==32) # 每个词在全图和区域级各计一次。
fig,ax=plt.subplots(figsize=(11,3.7));ids=np.arange(len(pa));ax.bar(ids-.17,pa,.34,label='layout A',color=BLUE);ax.bar(ids+.17,pb,.34,label='layout B',color=GOLD)
ax.set(xticks=np.arange(len(regions))*3+1,xticklabels=regions,xlabel='Region address (three A/B/C counts each)',ylabel='Count',title=f'Global + {gr} x {gc} regional counts: {len(pa)} dimensions')
for edge in np.arange(len(regions)+1)*3-.5:ax.axvline(edge,color='#ADBFD2',lw=.7)
ax.legend();fig.tight_layout();plt.show()
distance2=int(np.square(pa-pb).sum())
print('Squared distance: bag =',int(np.square(bag_a-bag_b).sum()),'; pyramid =',distance2)
RESULTS['pyramid']={'grid':[gr,gc],'dimensions':len(pa),'squared_distance':distance2,'vector_A':pa.tolist(),'vector_B':pb.tolist()}
print('Checks passed:',sum(CHECKS.values()),'/',len(CHECKS))

# 独立构造：只交换已提取描述的位置，词编号固定。
from collections import Counter
tokens=[('A',0.1),('A',0.2),('B',0.7),('B',0.8)]
swapped=[(word,1-y) for word,y in tokens]
def describe(rows):
    return dict(global_counts=dict(Counter(w for w,y in rows)),
                top=dict(Counter(w for w,y in rows if y<.5)),
                bottom=dict(Counter(w for w,y in rows if y>=.5)))
before,after=describe(tokens),describe(swapped)
print('原布局',before,'\n交换位置',after)
assert before['global_counts']==after['global_counts']
assert before['top']!=after['top']
print('HOG 配置维度：',(64//8-1)*(128//8-1)*2*2*9)
print('三层金字塔 K=3：',(1+4+16)*3)
assert abs(.16927+.19241-.36168)<1e-10